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SECTION  I 


INTRODUCTION 

Recent  interest  in  developing  high-speed  commercial  transports,  as  well  as  more  advanced 
reentry  vehicles,  has  lead  to  a  resurgence  of  research  into  hypersonic  aerodynamics.  Simulating 
this  flow  regime  requires  the  modeling  of  the  chemical  activity  that  characterizes  it.  Several  re¬ 
searchers  have  developed  successful  finite-rate  chemistry  codes  (see  Candler  and  MacCormack 
[Reference  1],  and  Walters  et  al.  [Reference  2]).  However,  finite-rate  chemistry  is  both  comph- 
cated  and  expensive.  In  addition  to  the  requirement  that  the  reaction  paths  be  known  (a  highly 
non-trivial  task  in  itself),  finite-rate  chemistry  requires  the  substitution  of  a  usually  large  num¬ 
ber  of  species  continuity  equations  in  place  of  the  global  mass  conservation  equation.  Also,  for 
near  equilibrium  flows,  the  finite-rate  equations  become  extremely  stiff.  Considering  the  diffi¬ 
culties  associated  with  finite-rate  chemistry,  flow  simulations  using  the  local  chemical  equilibri¬ 
um  assumption  become  very  attractive. 

In  recent  years,  the  numerical  solution  of  chemical  equilibrium  problems  has  been  investi¬ 
gated  by  several  researchers.  Liu  and  Vinokur  [Refo-ence  3]  give  an  excellent  review  of  algo¬ 
rithms  for  the  determination  of  properties  for  an  equilibrium  gas  mixture.  Meinqes  and  Morgan 
[Reference  4]  present  an  interesting  formulation  based  on  defining  all  the  necessary  reactions  in 
terms  of  element  variables. 

Many  researchers  have  develc^red  flow  solves  based  upon  curve  fits  for  the  equilibrium 
prq>erties  [Reference  3-9],  However,  it  needs  to  be  stressed  that  formulations  based  upon  curve 
fits,  while  competitive  in  terms  of  CPU  time,  are  limited  to  a  specific  mixture  composition  for 
each  tabulation  [Reference  10].  The  study  of  arbitrary  mixtures  requires  the  creation  of  new  tab¬ 
ulations  for  each  mixture  of  interest 

Other  researchers  have  explored  the  capabilities  of  a  flow  solver  based  on  the  local  clKmical 
equilibrium  behavior  of  a  given  gas  mixture.  Finite-rate  and  chemical  equilibrium  computations 
for  blunt-body  flow,  using  a  5-species  air  model,  were  compared  by  D6sid^,  Glinsky  and  Het- 
tena  [Reference  11].  Davy,  Lombard  and  Green  [Reference  12]  present  viscous  computations 
which  simulate  an  entry  into  the  atmosphere  of  Jupiter,  and  Wang  and  Chen  [Reference  13]  pres¬ 
ent  rocket  engine  nozzle  computations  using  hydrogen/air  mixtures.  However,  the  flow  solvers 
used  in  the  above  studies  are  incapable  of  handling  arbitrary  gas  mixtures,  and  the  codes  would 
have  to  be  rewritten  for  a  different  gas  mixture  of  interest. 

The  iq>proach  taken  in  this  study  involves  the  coupling  of  a  **Black  Box”  chemical  equilibii-  ■ 
um  solver  with  a  three-dimensional.  Thin  Layer  Navier-Stokes  flow  solver  that  has  bem  modi¬ 
fied  to  include  real  gas  effects.  The  Black  Box  computes  the  equilibrium  composition  and  tem- 
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perature  of  gas  mixtures  at  constant  density  and  internal  energy,  and  provides  the  flow  solver 
with  the  necessary  thermodynamic  and  transport  properties.  The  modifications  to  the  flow  solver 
are  minor,  although  they  include  the  implementation  of  a  newly  derived  approximate  Riemann 
solver  of  the  Roe  type.  A  key  advantage  of  this  “real  gas”  flow  solver  is  the  capability  to  handle 
arbitrary  mixtures  of  thermally  perfect  gases  in  local  chemical  equilibrium. 

The  subject  matter  of  the  present  study  is  divided  into  two  parts.  Part  one  concerns  the  de¬ 
velopment  and  testing  of  the  Black  Box  and  is  covered  by  Sections  n  through  V.  The  modifica¬ 
tion  of  a  perfect  gas  flow  solver  to  include  real  gas  effects  and  its  subsequent  application  com¬ 
prise  the  second  part  of  this  report,  which  is  covered  in  Sections  VI  through  Vin. 

The  mathematical  foundations  and  governing  equations  for  the  determination  of  the  equihb- 
rium  composition  and  thermodynamic  properties  of  a  homogeneous  mixture  of  thennally  perfect 
gases  are  developed  in  Section  n.  Two  different  solution  procedures  are  proposed,  and  a  third 
discussed  briefly.  In  Section  m,  practical  models  for  the  description  of  the  thermochemical 
properties  are  discussed.  Included  are  various  options  for  the  thermodynamic  model,  equilibrium 
constant,  chemistry  model,  and  transport  properties.  In  Section  IV,  numerical  procedures  for  the 
solution  of  the  governing  equations  developed  in  Section  n  are  given.  Results  from  the  numeri¬ 
cal  testing  of  the  Black  Box  chemical  equilitoium  solver  are  presented  in  Section  V.  Several 
chemistry  models  and  various  compositions  are  utilized,  and  efficiency  and  robustness  studies 
are  performed. 

The  governing  equations  for  three-dimensional,  viscous,  heat-conducting  flows  written  for 
a  time-dependent,  generalized  coordinate  system  are  presented  in  Section  VI.  In  Section  VII,  the 
numerical  methods  for  the  solution  of  the  gasdynamic  equations  are  presented.  Real  gas  formula¬ 
tions  are  developed  for  the  Steger- Wanning  flux-vector  splitting  scheme  and  split  flux  Jaco- 
bians,  as  well  as  a  newly  derived  approximate  Riemann  solver  of  the  Roe-type.  These  algo¬ 
rithms  are  cast  in  a  maimer  compatible  with  their  perfect  gas  count^arts.  Applications  of  the 
chemical  equilibrium  flow  solver  are  given  in  Section  Vm,  along  with  comparisons  with  perfect 
gas  solutions. 
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SECTION  II 


MATHEMATICAL  FORMULATION;  BLACK  BOX 
1.  CHEMICAL  EQUILIBRIUM 

In  order  to  classify  the  chemical  activity  that  characterizes  high-speed,  high-temperature 
flows,  an  order  of  magnitude  analysis  is  usually  made  of  the  reaction  time  versus  the  fluid  dy¬ 
namic  time.  Specifically,  considering  a  single  control  volume  in  the  flow  domain,  the  fluid  dy¬ 
namic  time  tpQ  can  be  defined  to  be  the  time  it  takes  for  a  fluid  particle  to  traverse  this  space. 
The  reaction  time  may  be  taken  as  the  time  required  to  partially  complete  a  chemical  reac¬ 
tion.  There  are  three  cases  that  occur. 


«  I^CB  * 

“tpj)*  ^CB  * 
»  ■'CB* 


(1) 


The  first  case  corresponds  to  frozen  flow,  in  which  the  reaction  has  no  time  to  occur.  The  second 
case  is  the  most  general  in  which  the  reaction  may  or  may  not  have  time  to  complete,  and  is  clas¬ 
sified  as  frnite-rate  chemistry.  For  the  last  case,  the  reaction  has  infinite  time  to  occur  and  thus 
will  reach  its  equilibrium  composition.  Physically,  frozen  and  equilibrium  chemistry  represent 
the  limiting  cases  for  the  real-life  situation:  frnite-rate  chemistry.  Furthermore,  if  the  last  case  is 
^plicable  at  each  volume  in  the  flowfreld  and  for  all  reactions  occurring  within  the  flow,  then 
the  flowfreld  can  be  assumed  to  be  in  local  chemical  equilibrium. 

In  the  following  development,  a  gas  mixture  composed  of  N  chemical  species,  in  which 
there  are  NE  elemental  species,  will  be  considered.  A  general  formulation  of  the  chemical  equa¬ 
tions  for  NR  reactions  involving  the  N  species  X|  can  be  written 

vij  Xi  +  V2J.  Xj  +  ...  +  vnj  Xf,  25 

v'lj  Xj  +  V2j  Xj  +  ...  +  Vn^  Xn,  r  =  1 . NR, 

where  vj^  are  the  stoichiometric  coefficients  of  species  i  in  reaction  r  for  the  reactants  and 
products  respectively. 

Considering  a  reactive  flow  involving  the  above  mentioned  gas  mixture,  the  governing  fluid 
dynamic  equations  can  be  derived  following  Vincenti  and  Kruger  [Reference  14].  In  particular, 
the  species  continuity  equations  written  in  integral  form  for  a  control  volume  V  bounded  by  a 
control  surface  S  are 
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Ttj  \  \  ^  I  -  Us)  ■  ndS  =  j  j  I  WidV  . 


1=1 . N,  (3) 


where  Qj  is  the  species  density,  Uj  is  the  species  velocity,  u*  is  the  surface  velocity  and  n  is  the 
unit  vector  normal  to  the  control  surface.  The  terms  in  the  above  equation  represent,  from  left  to 
right,  unsteady  contributions  integrated  over  the  volume,  inviscid  and  viscous  fluxes  integrated 
over  the  surface,  and  a  source  term  due  to  chemical  reactions  integrated  over  the  volume. 

TTie  source  term  w,  in  Equation  (3)  represents  the  rate  of  production  of  species  i  and  can  be 
written  as 

NR 

Wj  =  .ALi  ~  ^ir)  * 


i  =  1,...,N, 


where  is  the  molecular  mass  of  species  i  and  k^j. ,  k^j  are  the  forward  and  backward  reac¬ 
tion  rates  for  reaction  r,  respectively.  These  reaction  rates  are  related  by 

kfj.  =  k(,^cj-  ♦ 

where  the  equilibrium  constant  Kcj  is  introduced,  which  is  a  known  function  of  the  thermody¬ 
namic  state.  Substituting  the  preceding  relation  into  Equation  (4),  the  rate  of  production  can  be 
rewritten  in  the  following  form 


^n(€-nm 


i  =  1 . N. 


where  at  chemical  equilibrium  the  term  in  brackets  is  identically  zero  because  it  reduces  to  the 
Law  of  Mass  Action 


IT  —  1*1 

K-cj  -  “jf 


n(5:) 

1«1 

nfe) 


r  =  1,...,NR, 


which  is  valid  for  equilibrium  chemistry. 
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The  limiting  cases  of  equilibrium  and  frozen  chemistry  should  be  readily  apparent  in  Equa¬ 
tion  (6).  Frozen  chemistry  corresponds  to  the  backward  reaction  rate  going  to  zero,  which  results 
in  the  rate  of  production  terms  going  to  zero  and  dropping  from  the  species  continuity  equation. 
In  this  case  it  can  be  shown  that  there  will  be  no  change  in  the  species  mass  fractions  if  diffusion 
is  neglected.  On  the  other  hand,  equilibrium  chemistry  corresponds  to  the  backward  reaction 
rates  approaching  infinity.  In  this  case  the  bracketed  terms  in  Equation  (6)  go  to  zero  and  the  rate 
of  production  terms  will  reach  a  finite  limit  value  which  will  balance  out  the  unsteady,  inviscid 
and  viscous  contributions  to  the  species  continuity  equations.  Equation  (3). 

2.  EQUATIONS  OF  STATE 

A  very  important  assumption  in  the  following  development  is  that  the  effect  of  intermolecu- 
lar  forces  in  the  gas  mixture  is  negligible.  As  a  result,  the  individual  gases  comprising  the  mix¬ 
ture  will  behave  as  thermally  perfect  gases.  Moreover,  the  assumption  of  local  chemical  equilib¬ 
rium  implies  that  the  thermodynamic  state  at  any  given  time  and  any  poini  in  space  need  be 
defined  as  a  function  of  only  two  state  variables.  For  reasons  that  will  become  apparent  later, 
the  density  Q  and  temperature  T  were  chosen  as  the  two  fundamental  state  variables. 

For  a  mixture  composed  of  N  species,  the  total  density  q  and  the  mass  fractions  Yj  are  given 
by  the  relations 

o  =  Eei.  =  (8) 

i=l 

Using  the  standard  mixing  rule,  the  mixture  gas  constant  R  is  defined  from  the  known  species 
properties 

R=|;y,R,.  (9) 

i=l 

where  Rj  is  the  gas  constant  for  species  i.  Similarly,  the  total  internal  energy  of  a  mixture  can  be 
given  as  the  sum  of  the  species  contributions,  which  results  in  the  caloric  equation  of  state 

N  N 

e  =  2  Y^ei  =  ^  Yi  Cv,(x)  dt  +  h,^  ,  (10) 

i=l  i-i 

where  Cy,  is  the  species  specific  heat  at  constant  volume  and  hf  is  the  species  heat  of  fonnation 
at  reference  temperature  T^ . 

Since  the  mixture  is  composed  of  thermally  perfect  gases,  the  thermal  equation  of  state  will 
be  given  by  Dalton’s  Law 
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(11) 


N 

P  =  X  =  €>RT  , 

i=l 

where  the  total  mixture  pressure  is  the  sum  of  the  partial  pressures.  The  thermal  equation  of 
state  is  indirectly  related  to  the  caloric  equation  of  state  Equation  (10)  via  temperature.  Thus  for 
a  given  chemical  composition  the  temperature  must  be  obtained  from  either  the  caloric  or  the 
thermal  equation  of  state  depending  on  whether  internal  energy  or  pressure  is  known. 


The  mixture  specific  heats  can  be  derived  as  functions  of  density  and  temperature  deriva¬ 
tives  of  the  mass  fractions.  The  specific  heat  at  constant  volume  is  defined  to  be 


=  =  V  Y  / +  'Vi- ( 

~  laT/v  2-^4aT)  ■‘‘Z.MaT)  • 

i=l  '  'Q  i=l  '  'Q 


(12) 


Recognizing  that  the  first  partial  is  nothing  more  than  the  species  specific  heat  at  constant  vol¬ 
ume,  the  relation  will  reduce  to 


Cv  Cy 


where  the  frozen  specific  heat  at  constant  volume 

N 

Cv  =  ^  Y^Cy.  , 
i=l 

has  been  introduced.  The  mixture  specific  heat  at  constant  pressure  can  be  defined  as 


c 


p  “ 


(13) 


(14) 


(15) 


where  the  mixture  enthalpy  h  =  e  +  RT  has  been  introduced,  which  can  be  given  as  the  sum  of 
the  species  contributions 


h  = 


using  the  standard  mixture  rule.  The  enthalpy  derivatives  are  then  obtained  as 


/  \  N  /av-  \ 

lih]  -  Vhl£Iil 


and 


(16) 


(17) 


(18) 
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where  the  frozen  specific  heat  at  constant  pressure  has  been  introduced 


Cp 

i=l 

The  partial  derivative  of  density  with  respect  to  temperature  at  constant  pressure  can  be  obtained 
from  the  differential  form  of  the  thermal  equation  of  state,  Equation  (11),  as 


R  +  T 


=  _  £  i=i  ^  ^ 

[dT  I  T  N  /  V 


(20) 


Finally,  substitution  of  Equations  (17-20)  into  Equation  (15)  gives  the  mixture  specific  heat  at 
constant  pressure 


Cp  —  Cv  "h 


(21) 


3.  SPEED  OF  SOUND 


A  key  thermodynamic  prq)erty  for  high-speed  flows  is  the  speed  of  sound,  which  is  defined 
as 


a 


2 


(22) 


For  a  thermodynamic  system  in  chemical  equilibrium,  the  combined  First  arxl  SecoiKl  Laws  of 
Theimodynamics  can  be  written  as 


Tds  =  de 


(23) 


Using  the  above  relation,  the  equation  for  the  speed  of  sound  Equation  (22)  can  be  rewritten  in 
the  following  form 


Using  the  differential  forms  of  the  caloric  and  thermal  equations  of  state.  Equations  (10)  and 
(11)  respectively,  and  Equation  (23),  the  pressure  partial  derivatives  are  found  to  be 
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Substituting  Equations  (25)  and  (26)  into  Equation  (24),  the  equation  for  the  speed  of  sound  can 
be  put  into  a  more  familiar  form 

=  r  I ,  (27) 

where  the  isentropic  index  F  is  defined  to  be 

r  =  (28) 

In  the  above  equation,  the  ratio  of  the  partial  derivative  of  enthalpy  with  respect  to  tempo-ature 
at  constant  density  to  the  partial  derivative  of  energy  with  respect  to  temperature  at  constant  den¬ 
sity  has  been  introduced 


Although  the  above  equation  for  the  speed  of  souikI  is  seemingly  simple,  it  should  be  readily  ap¬ 
parent  that  the  isentropic  index  is  a  complex  fimction  of  density  and  temperature. 

In  general,  there  are  at  least  four  different  ‘‘gammas”  that  can  be  defined  for  an  equilibrium 
gas:  the  two  mentioned  already,  the  ratio  of  specific  heats  y  =  Cp/cv,  and  the  ratio  of  frozen 
specific  heats  y  =  Cp/cy.  Although  generally  their  values  are  different,  all  four  of  these  “gam¬ 
mas”  will  reduce  to  the  same  value  for  firozen  flows,  of  which  pofect  gas  is  a  subset  More  de¬ 
tails  on  the  derivation  of  the  equilibrium  sound  speed  are  given  in  Appendix  A. 
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4.  SOLUTION  TECHNIQUES 


As  stated  previously,  chemical  equilibrium  is  achieved  when  the  reaction  rates  go  to  infiraty. 
Thus  the  species  rate  of  production  terms  reach  a  limiting  value  and  do  so  while  keeping  two 
state  properties  constant.  Examination  of  the  species  continuity  equations.  Equation  (3),  shows 
that  changes  in  these  production  terms  occur  at  constant  specific  volume  and  thus  constant  densi¬ 
ty.  From  the  combined  First  and  Second  Laws  of  Thermodynamics,  Equation  (23),  it  can  be  seen 
that  for  a  reversible  process  occurring  at  constant  density,  the  internal  energy  will  also  remain 
constant.  Therefore,  evaluation  of  the  equilibrium  system  will  occur  at  constant  density  and  in¬ 
ternal  energy. 

Two  different  formulations  for  the  solution  of  the  chemical  equilibrium  state  at  constant 
density  and  internal  energy  have  been  developed  and  will  be  discussed  in  detail  in  the  following 
subsections.  The  first  is  the  Mass  Constraint  Technique,  which  uses  a  modified  form  of  the  laws 
of  mass  action,  plus  elemental  mass  constraints  and  an  energy  equation.  The  resulting  system  of 
equations  can  be  reduced  through  the  use  of  the  concept  of  degree  of  advancement  for  a  chemi¬ 
cal  reaction.  This  reduction  is  the  foundation  of  the  second  formulation,  which  will  be  named  the 
Degree  of  Advancement  Technique.  Additionally,  a  third  formulation,  based  on  the  work  of 
Meintjes  and  Morgan  [Reference  4],  was  attempted.  However,  preliminary  tests  indicated  that 
this  formulation  was  not  suitable  for  calculations  where  temperature  is  not  constant. 

In  addition  to  the  previous  formulations,  the  first  method  has  been  modified  in  order  to 
solve  chemical  equilibrium  problems  at  constant  pressure  and  temperature,  and  constant  pressure 
and  density.  The  importance  of  these  modifications  will  become  apparent  later,  when  the  bound¬ 
ary  conditions  for  the  flow  solver  will  be  described. 

a.  MASS  CONSTRAINT  TECHNIQUE 

The  chemical  equilibrium  state  can  be  determined  through  the  solution  of  an  algebraic  set  of 
N  -I-  1  equations,  comprising  N-NE  laws  of  mass  action,  NE  elemental  species  mass  constraints 
and  an  energy  equation,  for  N  -I-  1  unknowns,  N  species  densities  and  temperature.  In  this  devel¬ 
opment,  the  governing  equations  are  written  in  a  form  compatible  with  finite-rate  chemistry. 

As  stated  previously,  the  limiting  case  of  chemical  equilibrium  corresponds  to  backward 
reaction  rates  going  to  infinity,  and  the  bracketed  terms  in  Equation  (6)  going  to  zero  since  they 
reduce  to  the  laws  of  mass  action.  In  order  to  avoid  the  difficulties  associated  with  the  deter¬ 
mination  of  the  limit  values,  the  backward  reaction  rates  are  assumed  to  have  a  firate  value,  uni¬ 
ty,  and  the  rates  of  productions  are  set  equal  to  zero.  This  gives  a  modified  form  of  the  laws  of 
mass  action  that  remains  compatible  to  the  finite-rate  form  expressed  in  Equation  (6)  while  still 
satisfying  the  chemical  equilibrium  state.  These  modified  laws  of  mass  action  are  written  as 
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=  0  , 


(30) 


NR 

Wj  (Qj,T)  =  -  vij) 

r=l 


where  the  index  i  now  covers  only  the  nonelemental  species  instead  of  the  entire  range. 

The  elemental  mass  conservation  equations  represent  the  fact  that,  in  the  absence  of  nuclear 
reactions,  mass  is  neither  created  nor  destroyed,  and  are  written  as 


'N-NE  +  i 


. NE- 

j=l  ^  U  =  l  J  Jt=n 


where  a^  is  the  number  of  particles  of  element  i  present  in  species  j.  The  second  term  on  the  left 

hand  side  represents  the  initial  values  at  the  beginning  of  the  computation.  The  NE  elemental 
species  need  not  be  true  elements  themselves,  but  they  must  contain  the  element  that  they  repre¬ 
sent.  An  example  would  be  the  use  of  carbon  monoxide  CO  as  the  elemental  species  for  either 
the  element  carbon  C  or  the  element  oxygen  O. 


In  order  to  complete  the  system,  an  equation  relating  the  unknown  temperature  to  the  known 
state  properties  is  needed.  This  is  accomplished  with  the  caloric  equation  of  state.  Equation  (10), 
written  in  terms  of  the  total  energy  per  unit  volume 


d(Qj.T)  =  Qeo  -  ^Qj  Cy/t)dt  -I-  hf. 


where  density  q,  total  energy  per  unit  volume  peQ,  and  velocity  Uq  are  known  quantities  that  are 
readily  available  from  any  flow  solver  written  in  conservation  form.  The  species  specific  heats  at 
constant  volume  Cy^  are  functions  of  temperature  and  are  provided  by  the  thermodynamic  mod¬ 
els,  which  win  be  discussed  in  Section  m. 


b.  DEGREE  OF  ADVANCEMENT  TECHNIQUE 

The  use  of  the  concq)t  of  degree  of  advancement  for  a  chemical  reaction  provides  an  alter¬ 
native  formulation  to  the  Mass  Constraint  Technique.  This  alternative  formulation,  titled  the  De¬ 
gree  of  Advancement  Technique,  can  provide  a  significant  advantage  over  the  previous  formula¬ 
tion.  The  advantage  comes  from  reducing  the  system  of  N+1  equations  to  as  little  as  N-NE+1 
equations  by  substituting  degrees  of  advancements  for  the  species  densities. 

For  a  single  reaction  of  the  form  given  by  Equation  (2),  the  ratio  of  the  infinitesimal 
changes  in  the  number  of  moles  for  each  species  is  proportional  to  the  ratio  of  the  differences 
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between  the  product  and  reactant  stoichiometric  coefficients  of  those  species.  This  relationship  is 
written  as 


dJTi  :  dJTj  :  ...  :  =  (v’,’  -  vl) :  (vj’  -  vj) : ...  :  (vn  -  Vn)  ,  (33) 

where  JTj  is  the  mole  number  for  the  species  j.  An  example  of  this  relationship  can  be  given  for 
the  simple  combustion  reaction 

2H2  +  02^  2H2O  ,  (34) 

where  Equation  (33)  will  read 


:  dJTo^ :  =  (-  2) ;  (-  1) :  (+  2)  .  (35) 

Introducing  an  infinitesimal  constant  of  proportionality  and  recognizing  that  dN'j  =  dYj/.4t)j, 
where  Yj  =  Qj/q  is  the  species  mass  fraction,  these  ratios  can  be  rewritten  as 

dY-  dYx, 

’x  =  -  =  ^  .  (36) 

,4bj(Vi  Vj)  ./H32(V2  V2)  —  Vjj) 

where  |  is  the  degree  of  advancement  for  the  reaction.  Extending  this  ^proach  to  multiple, 
coupled  reactions,  the  degree  of  advancement  for  reaction  r  can  be  defined  by  the  differential 
relation 


=  .4bj(v;;  -  vj^)  ,  j  =  1,...,N,  r  =  1 . NR  .  (37) 

Noting  that  density  is  kept  constant,  the  system  of  NR  equations  for  a  species  j  can  be  solved  to 
give  the  species  density  as  a  function  of  the  NR  degrees  of  advancement 


NR 

Oj  “  (Qpt=0  Q  ~  Yj,r)^r  ♦  j  ~  1,...»N,  (38) 

r=l 

where  the  initial  condition  (|r)t»o  =  (^  been  applied. 

Using  the  above  transformation  relation.  Equation  (38),  the  laws  of  mass  action  as  well  as  the 
energy  equation  can  be  rewritten  in  terms  of  the  degrees  of  advancement  Furthermore,  substitut¬ 
ing  Equation  (38)  into  the  mass  constraint  equations  Equation  (31),  as  follows 
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results,  after  some  manipulation,  into  the  expression 

NR 

r=l 

It  is  easily  verified  that  the  bracketed  terms  go  to  zero  identically  since  reaction  r  must  be  stoi- 
chiometrically  balanced  with  respect  to  element  i.  To  illustrate  this,  consider  again  the  example 
of  the  simple  combustion  reaction  consisting  of  the  three  species  H2, 02  and  H2O.  The  elemental 
species  are  selected  as  H2  and  O2  and  represent  the  two  elements,  atomic  hydrogen  H  and  oxy¬ 
gen  O,  respectively.  Expanding  the  bracketed  terms  for  each  element  gives 

2(0  -2)  +  0(0  -  1)  +  2(2  -  0)  s  0,  Element  =  H,  (41) 

0(0  -2)  +  2(0  -  1)  +  1(2  -  0)  s  0,  Element  =  O  .  (42) 

Thus  the  mass  constraints  are  implicitly  satisfied  and  are  eliminated  firom  the  governing  equa¬ 
tions  because  the  degrees  of  advancement  will  stoichiometiically  balance  the  individual  reac¬ 
tions.  The  resulting  system  of  equations  consists  of  NR  laws  of  mass  action  and  an  energy  equa¬ 
tion,  with  the  unknowns  being  the  degrees  of  advancement  and  temperature.  The  maximum 
advantage  over  the  Mass  Constraint  Technique  is  obtained  when  the  number  of  reactions  is  equal 
to  the  minimum  number  of  necessary  reactions,  NR  =  N  -  NE  . 

c.  SYSTEMS  AS  FUNCTIONS  OF  ALTERNATE  STATE  PROPERTIES 

In  the  previous  sections  the  caloric  equation  of  state  was  employed  to  find  the  chemical 
composition  of  the  gas  mixture.  This  choice  works  well  for  control  volumes  in  the  flow  field, 
since  the  two  known  state  variables  are  density  and  internal  energy.  However,  either  or  both  of 
these  state  variables  may  not  be  available  at  the  boundaries  if  certain  boundary  conditions  are  to 
be  enforced.  For  example,  the  known  state  variables  for  specified  pressure  and  temperature  gra¬ 
dient  boundary  conditions  are  pressure  and  temperature.  Another  example  is  subsonic  character¬ 
istic  variable  boundary  conditions  whereby  pressure  and  density  are  known.  Consequently,  the 
determination  of  the  equilibrium  composition  at  boundary  points  has  to  be  handled  separately. 
When  pressure  is  a  known  property,  as  in  both  of  the  above  mentioned  examples,  the  thermal 
equation  of  state.  Equation  (11),  can  be  employed  instead  of  the  caloric  equation  of  state.  The 
other  equations  in  the  system  remain  unchanged.  However,  while  the  unknowns  are  still  species 
densities  and  temperature  for  the  first  example,  (p  —  Q  case),  they  become  the  species  densities 
and  total  density  for  the  second  example,  (p  -  T  case). 


N 

E  -  Vj’^)  I,  =  0,  1=1 . NE  .  (40) 

i=i 
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5.  THERMODYNAMIC  PROPERTIES 


The  next  step  following  the  determination  of  the  equilibrium  composition  is  the  ev  luation 
of  any  thermodynamic  properties  of  interest,  which  include  pressure,  speed  of  sound,  and  isen- 
tropic  index  F  (commonly  used  in  flux-split  numerical  schemes). 

Evaluation  of  some  of  these  properties  requires  knowledge  of  partial  derivatives  of  the  mass 
fractions  with  respect  to  density  and  temperature,  as  can  be  seen  from  the  expression  for  the 
speed  of  sound.  Equations  (27-29).  These  doivatives  can  be  expressed  as  functions  of  the  partial 
derivatives  of  the  species  densities,  and  are  given  as 


The  derivatives  of  the  species  densities  can  be  obtained  using  the  Imphcit  Function  Theo¬ 
rem.  Specifically,  the  first  set  of  N  equations  in  the  M-ss  i.,onstfamt  Technique  implicitly  define 
N  functions  of  density  and  temperature 

Wi[Qi(e,T),e/Q,T),...,eN<Q,T)]  =  o  ,  i  =  i,...,n  ,  (45) 

where  the  species  densities  are  the  N  functions,  Pj  =  Qj(Q,  T).  The  Implicit  Function  Theorem 
states  that  the  derivatives  of  the  N  functions  with  respect  to  density  and  temperature  are  given  by 
the  solution  of  dwj  =  0  provided  that  the  partial  derivatives  of  w^  exist  and  are  continuous,  and 
the  determinant  of  matrix  A  =  dWj/dQj  is  not  zero.  Thus  the  derivatives  of  the  species  densities 

with  respect  to  density  and  temperature  are  obtained  through  the  solution  of  the  following  sys¬ 
tems 


where  B  =  d'f/  Jdq  and  C  =  dy/JdT. 

This  procedure  can  be  similarly  applied  to  the  E>egiee  of  Advancement  Technique,  where 
the  functions  are  now  N  -  NE  degrees  of  advancements,  |r  =  lr(Q.  T).  The  linear  systems  are 
now 
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(48) 


where  matrix  A  =  3wj/ and  the  right-hand  sides  of  both  systems  remain  the  same.  Using 
the  differential  form  of  tire  transformation  relation.  Equation  (37),  the  partial  derivatives  of  the 
mass  fractions  can  be  given  in  terms  of  the  partial  derivatives  of  the  degrees  of  advancement,  as 
follows 


(50) 

(51) 
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SECTION  m 


PRACTICAL  MODELS 

The  mathematical  foundation  of  local  chemical  equihbrium  has  been  given  in  the  previous 
section  for  an  arbitrary  mixture  of  thermally  perfect  gases.  No  restrictive  assumptions  have  been 
placed  on  the  functional  form  of  the  species  internal  energy,  with  the  exception  that  each  compo¬ 
nent  is  considered  to  behave  as  a  thermally  perfect  gas.  Similarly,  the  equilibrium  constant,  the 
chemical  system  and  the  other  related  properties  have  been  given  in  general  form.  Therefore, 
virtually  any  practical  model  can  be  utilized  in  conjunction  with  the  present  formulation.  The 
specific  models  developed  and  used  in  the  “Black  Box”  for  the  species  internal  energy,  the  equi¬ 
librium  constant,  the  chemical  system  and  the  transport  properties  are  given  in  the  following. 

1.  INTERNAL  ENERGY 

There  are  two  thermodynamic  models  available  for  use  in  the  Black  Box.  The  first  is  an 
equilibrium  statistical  mechanics  model  and  will  be  referred  to  as  the  Vibrational  Model.  The 
second  model  is  based  on  polynomial  curvefits  and  will  be  termed  the  Curvefit  Model. 

a.  VIBRAnONAL  MODEL 

From  statistical  mechanics,  the  internal  energy  can  be  taken  as  the  sum  of  the  translational, 
rotational,  vibratioiuil  and  electronic  contributions  [Reference  14].  Neglecting  the  electronic  con¬ 
tributions  and  considering  the  translational  and  rotational  contributions  at  their  fully  excited  val¬ 
ues  the  equation  for  the  species  internal  energy  will  read 


where  the  vibrational  contributions  are  modeled  as  simple  harmonic  oscillators.  In  the  above, 
are  the  characteristic  temperatures  for  vibration  arxl  the  summation  represents  the  fact  that 
for  polyatomic  molecules  there  will  be  more  than  one  vibrational  temperature.  The  number  of 
vibrational  temperatures  for  species  i  is  defined  to  be  NVT{.  In  Equation  (S2),  the  first  term  rq)- 
resents  the  translational  and  rotational  contributions.  The  parameter  nt  will  have  a  value  of  3/2 
for  a  monatomic  gas,  5/2  for  a  diatomic  or  linear  polyatomic  gas  and  6/2  for  a  nonlinear  poly¬ 
atomic  gas.  The  characteristic  vibrational  temperatures  are  obtained  fi-om  the  JANAF  tables 
[Reference  22]. 
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b.  CURVEFTT  MODEL 


The  Curvefit  Model  uses  fourth-order  polynomials  to  curvefit  the  species  specific  heat  at 
constant  volume 

Cvj  ~  8^14  "b  ^24  %4  ^44  ^54  »  (53) 

where  the  internal  energy  is  given  by 

®i  =  *^64+  I  Cvi(t)dt  ,  i  =  (54) 

There  are  several  curvefit  formulations  available  in  the  literature.  McBride  et  al  [Reference  15]. 
give  curvefits  for  molecules  involving  the  first  18  elements,  which  are  iq)plicable  up  to  6000  K. 
This  range  was  extended  to  15,000  K  for  carbon,  hydrogen,  nitrogen  and  oxygen  systems  by 
Esch  et  al  [Reference  16].  Gupta  et  al  [Reference  17].  provide  improved  curvefits,  ^plicable  up 
to  30,000  K  for  an  1 1-species  air  model.  In  the  Black  Box,  the  first  set  of  curvefits  is  used  for 
argon  systems,  the  second  set  for  carbon  and  hydrogen  systems,  and  the  latter  curvefits  are  uti¬ 
lized  for  the  eleven  species  air  model  they  were  intended  for. 

2.  EQUILIBRIUM  CONSTANT 

The  equilibrium  constants  define  the  equilibrium  composition  of  the  mixture  and  are  fimc- 
tions  of  temperature.  The  Black  Box  incorporates  two  models  for  their  determination.  In  the  first 
one,  experimental  values  are  curvefitted  employing  an  Amiehius-type  functional  form.  The  se¬ 
cond  model  derives  the  equilibrium  constant  from  the  species  Gibb’s  free  energy,  following  the 
example  of  liepmatm  and  Roshko  [Reference  18]. 

a.  CURVEFIT  Kc 

The  first  model  is  based  on  curvefits  compatible  with  Arrenhius  type  chemistry 

Kcj  =  Cc'n«e"®'^  ,  r  =  1 . NR  ,  (55) 

where  the  constants  can  be  found  in  Vincent!  and  Kruger  [Reference  14]  for  a  7-species  air  mod¬ 
el  and  Kang  and  Dunn  [Reference  19]  for  an  11-species  air  model.  The  major  advantage  of  the 
Curvefit  Kc  is  the  use  of  accurate  empirical  data  for  the  determination  of  the  equilibrium 
constant  However,  this  leads  to  a  disadvantage  in  tfiat  the  empirical  data  may  not  be  available. 
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b.  CONSISTENT  Kc 


The  second  model  is  based  on  the  minimization  of  Gibb’s  free  energy  at  constant  pressure 
and  temperature.  Thus  the  equilibrium  constant  can  be  obtained  indirectly  from  the  thermody¬ 
namic  model  and  therefore  is  not  limited  to  the  availability  of  experimental  data.  This  model, 
termed  the  Consistent  Model,  is  very  general  and  depends  only  on  the  availability  of  certain 
reference  values  which  are  much  easier  to  obtain  than  the  curvefits  mentioned  previously.  How¬ 
ever,  any  errors  introduced  by  assumptions  and  simplifications  made  in  the  thermodynamic  mod¬ 
el  are  inherited  by  the  equilibrium  constant. 

Gibbs’s  free  energy  for  a  species  can  be  written  as 


gi  =  hi  -  Tsj  , 
T 


=  [  c„(t)dt  +  h,,  -  T  [  +  R,Tln(|i)  -  Ts,,  , 


(56) 


where  the  equation  for  species  entropy  and  enthalpy  have  been  utilized.  In  the  above,  s,  and  h, 
are  the  reference  entropy  and  enthalpy  taken  at  the  reference  temperature  T^f  and  pressure  Pf . 
Using  the  standard  mixture  rule,  the  total  Gibb’s  free  energy  for  the  mixture  is  then 


Gibb’s  free  energy  can  be  written  as  a  function  of  pressure,  temperature  and  the  degrees  of  ad¬ 
vancement,  g  =  g(  p,T,|r).  Chemical  equilibrium  occurs  at  the  point  where  the  Gibb’s  free 
energy  is  a  minimum  under  constant  pressure  and  temperature.  Minimizing  the  above  equation 
for  reaction  r  at  constant  pressure  and  temperature,  dg/d^r  =  0,  utilizing  the  differential  form  of 
the  transformation  relation  Equation  (37),  and  taking  the  exponential  of  both  sides,  results  in  the 
expression 


where 


“•  =  gi  -  RiTln 


[  Cp,(t)dt-T[  +  h,  _  Ts,_  . 


(59) 
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The  previous  formula  is  the  Law  of  Mass  Action  written  in  terms  of  the  partial  pressures.  Ex¬ 
pressing  the  constant  in  terms  of  the  concentrations  and  assuming  that  the  reference  pressure  is 
the  same  for  all  species  (Pf^  =  p^^f),  the  following  relation  for  the  equilibrium  constant  is  ob¬ 
tained 


K 


CJ 


=  exp 

\rt/  ' 


N 


1=1 


(60) 


where 


a.(T)  =  5^  .  (61) 

The  term  o),  given  in  Equation  is  determined  by  the  thermodynamic  model  employed.  In  the  fi¬ 
nal  form,  the  term  Qg  will  be  given  for  the  Curvefit  Model  as  follows 

a.(T)  =  a,,(l  -  InT)  -  ^  ,  (62) 

and  for  the  Vibrational  Model  as 


Q,(T)  =  l-l-  n,  +  ln(l  -  -  (1  +  n,)lnT  +  hJ/R.T  -  s°/R,  .  (63) 

The  last  two  terras  of  the  above  relations  are  identical,  because  the  constants  a^^  and  a7^  read 

a64  =  h?/R.  ,  a7^  =  s?/R.  ,  (64) 

where  h^  and  s?  are  the  species  enthalpy  and  entrq}y  at  a  reference  temperature  of  absolute  zero. 


3.  CHEMISTRY  MODELS 


As  stated  in  Section  I,  a  main  requirement  of  the  solver  is  the  capability  to  haixUe  arbitrary 
mixtures.  The  incorporation  of  additional  models  to  the  Black  Box  is  a  relatively  painless  pro¬ 
cess,  provided  that  a  few  key  reference  values  (i.e.  heats  of  formation)  can  be  supplied  for  the 
mixture  components.  The  Black  Box  currently  has  available  12  diffoent  chemistry  models. 
These  include  6  air  models  (ranging  firom  a  simple  5-species  to  a  relatively  complex  17-spe¬ 
cies),  3  combustion  models  arid  an  argon  plasma  model.  In  order  to  provide  the  flow  solver  with 
the  c^ability  to  revert  to  a  perfect  gas  solver,  a  perfect  gas  model  is  also  included,  with  the  spe¬ 
cific  gas  being  selected  at  input  from  a  set  comprising  air,  Ar,  CO2,  CO,  N2, 02>  H2  and  steam. 
Debugging  and  testing  of  the  code  was  accomplished  using  an  ideally  dissociating  oxygen  model 

OjS?20,  (65) 

where  the  constants  for  the  dissociation  of  oxygen  were  obtained  from  V^ncenti  and  Kruger[Ref- 
erence  14].  In  the  following,  the  latter  model  will  be  referred  to  as  Oxygen  Model. 
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a.  AIR  MODELS 


The  first  air  model  is  composed  of  five  species,  N2, 02,  NO,  N  and  O,  where  three  indepen¬ 
dent  reactions  are  used 


Nj  25  2N  , 

(66) 

O2  ^  20  , 

(67) 

NO  25  N  +  0  . 

(68) 

Atomic  nitrogen  and  oxygen  are  the  two  elemental  species  and  N2, 02,  NO  are  the  non-elemen- 
tal  species.  Two  mass  constraints  are  used  to  preserve  the  total  mass  of  atomic  oxygen  aixl  nitro¬ 
gen.  This  model  will  be  referred  to  as  S-Species  Air  Model.  Unless  otherwise  noted,  the  air 
models  will  be  lefened  to  by  the  number  of  species  present 

The  second  model  adds  two  more  species,  NO'*'  and  e~,  to  the  model  listed  above,  for  a 
total  of  seven.  The  additional  independent  reaction 

NOisNO'*' +  e-  ,  (69) 

is  needed  to  describe  the  ionization  of  NO.  Conservation  of  the  total  number  of  electrons  is  pro¬ 
vided  as  an  additional  mass  constraint  and  NO'*'  will  act  as  a  non-elemental  species. 

The  tfiird  model  utilizes  nine  species,  adding  N***  and  O'*'  to  the  seven  above.  Two  more 
indq>endent  reactions  are  needed,  describing  the  ionization  of  atomic  oxygen  and  nitrogen 

N2?N'''+e",  (70) 

0  25  0'''+e",  (71) 

and  bringu^  the  total  number  of  reactions  to  six.  No  additional  mass  constraints  are  necessary, 
due  to  the  fact  that  the  two  additional  species  are  non-elemental  species. 

Two  more  iixlqpendent  reactions  are  needed  for  the  fourth  model 

N2  N+  +  e"  ,  (72) 

02  25  02++e-,  (73) 

which  comprises  eleven  species,  with  the  addition  of  N^  and  to  die  previous  nine.  Again, 
there  are  no  new  mass  constraints,  because  the  additional  species  are  not  elements. 

The  fifth  model  is  obtained  by  adding  the  base  elemoit  argon  to  the  previous  model,  for  a 
total  of  thirteen  species.  The  two  additional  species  considered  are  Ar  and  Ar^.  The  ionization 
of  argon  is  described  by  the  additional  reaction 
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At  2^  Ar"^  +  e  , 


(74) 


and  Ax'*'  acts  as  a  nonelemental  species.  An  additional  mass  constraint  for  the  preservation  of 
argon  is  also  required. 

The  sixth  and  most  complex  air  model  employs  seventeen  species,  the  thirteen  above  plus  C, 
C  CO  and  CO2.  Three  more  indqrendent  reactions  are  needed 

2CO2  2CO  +  O2  ,  (75) 

2Ca^  20  +  02,  (”76) 

C2¥C-^+e-,  (77) 

along  with  die  previous  nine.  The  inclusion  of  carbon  requires  the  addition  of  anoth^  mass 
constraint  conserving  this  element,  and  the  additional  nonelemental  species  are  C  ,  CO  and 
CO2.  The  17-Species  Air  Model  employs  consistent  equilibrium  constants  only,  since  no  curve- 
fits  for  the  last  three  equations  have  been  fouiul. 

b.  COMBUSTION  MODELS 

Combustion  mixtures  provide  good  examples  of  the  difficulties  associated  with  the  use  of 
curvefits  for  the  determination  of  equilibrium  ccmipositions,  since  th^e  are  nearly  infinite  possi¬ 
bilities  for  the  mixture  ratio.  Utilizing  curvefits  would  require  a  matching  number  of  tabulations. 
The  only  viable  alternatives  are  either  the  mediod  presented  in  this  study  or  finite-rate  chemistry 
investigations  with  all  their  inherent  complexities. 

The  first  combustion  model  is  the  simple  two-reaction  hydrogen  combustion  in  air  pres¬ 
ented  in  Rogers  and  Chinitz  [Reference  20],  and  will  be  refenred  to  as  the  Hydrogen-Air  Com¬ 
bustion  Model.  A  total  of  five  species  are  utilized,  where  the  nonelemental  species  are  H2  and 
O2.  The  two  independent  reactions  are 


H2  +  O2  20H  , 

(78) 

20H  +  O2  25  2H2O  , 

(79) 

where  nitrogen  is  considered  to  be  inert,  aixl  three  mass  ccxistraints  provide  for  the  conservation 
of  nitrogen,  oxygen  and  hydrogen.  N2,  OH,  and  H2O  were  selected  as  the  elem^ital  species, 
where  the  latter  two  were  chosen  in  order  to  show  how  arbitrary  die  differentiation  betwem  ele¬ 
mental  and  nonelemental  species  can  be. 

The  second  combustion  model  corresponds  to  the  combustion  of  hydrogen  in  oxygen  arxl 
will  be  referred  to  as  Hydrogen-Oxygen  Combustion  Model.  The  six  species  O2,  H2,  OH,  H2O, 
O  and  H  comprise  the  mixture  and  the  four  required  reactions  are 
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02^20  , 

(80) 

Hj  25  2H  , 

(81) 

OH  25  0  +  H  , 

00 

s> 

HjO  25  2H  +  0  , 

(83) 

Here,  two  mass  constraints  are  used  for  the  preservation  of  oxygen  and  hydrogen,  and  the  ele¬ 
mental  species  are  O  and  H. 

The  last  model  is  taken  from  Mein^  and  Morgan  [Reference  4],  and  is  composed  of  the 
above  six  species  plus  CC)2i  CO  and  N2.  One  additional  independent  reaction  is  required 

CO2  CO  +  O  .  (84) 

Conservation  of  hydrogen,  oxygen,  carbon  and  nitrogen  is  ensured  dirough  the  use  of  four  mass 
constraints,  where  the  elemental  species  are  H,  O,  CO  and  die  inert  N2.  This  model  will  be  re¬ 
ferred  to  as  the  Hydrocarbon  Combustion  Model. 

c.  ARGON  PLASMA  MODEL 

The  Argon  Model  provides  a  simple  example  of  high  energy  plasmas.  It  employs  three  spe¬ 
cies  At'*',  Ar  and  e~,  and  one  reaction  is  required 

Ar  2*  Ar'*'  +  e”  ,  (85) 

describing  the  ionization  of  argon.  Two  mass  constraints  are  employed  for  the  conservation  of 
electrons  and  argon,  and  Ar'*'  will  act  as  the  ncmetemental  species.  Constants  for  the  equilibrium 
constant  are  obtained  from  Glass  and  Takano  [Reference  21]. 

Implementation  of  all  the  above  models  in  the  Black  Box  is  accomplished  dirough  the  use  of 
a  chemical  model  database.  For  each  of  die  species  listed  above,  dns  database  contains  all  the 
various  chemical  and  thennodynamic  properties  necessary  for  the  solution  of  the  governing 
equations.  Moreover,  the  stoichiometric  coefficients  and  curvefit  constants  for  various  reactions 
are  given  in  terms  of  all  the  species  available  in  the  database.  The  chonistry  model  of  interest  is 
then  constructed  from  this  database  by  specifying  which  species  are  presmit,  which  reactimis  to 
use,  etc.  The  presence  of  the  database  makes  the  implementaticm  of  additional  models  a  very 
sirr^le  process. 

4.  TRANSPORT  PROPERITES 

After  the  equilibrium  ccxnposition  has  been  obtained  by  means  of  the  techniques  discussed 
in  Section  L  the  mixture  transport  properties  can  be  easily  evaluated  by  following  two  sequential 
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steps.  The  first  is  the  determination  of  the  properties  for  the  individual  species,  whereas  the  se¬ 
cond  involves  the  application  of  a  mixing  rule  to  obtain  the  global  values. 


a.  SPECIES  VISCOSITY 


Two  models  for  species  viscosity  are  implemented  in  the  Black  Box.  The  first  model  for 
species  viscosity  is  a  straightforward  extension  of  Sutherland’s  Law  to  multicomponent  mixtures 
and  is  referred  to  as  the  Sutherland  Model.  The  species  viscosity  is  given  by  the  relation 


= 


t  +  Gm  ’ 


(86) 


where  the  constant  and  are  empirically  determined. 

The  second  model  is  based  on  the  curvefit  tabulations  by  Gupta  et  al.  [Reference  17],  which 
are  valid  for  air  chemistry  up  to  eleven  species  and  30,000  K.  The  equation  for  species  viscosity 
will  read 


(X-  =  ,  (87) 

where  B^  and  are  the  curvefit  coefficients.  This  model  will  be  referred  to  as  Gupta 
Model.  While  more  accurate  at  higher  temperatures  than  Sutherland’s  Law,  the  Gupta  Model 
does  not  fair  well  below  1(XX)  K.  Therefore  for  temperatures  below  this  limit,  the  Sutherland 
Model  will  be  used. 


b.  SPECIES  THERMAL  CONDUCTIVITY 


Similarly  to  what  has  been  seen  for  viscosity,  two  models  are  employed  for  the  species  ther¬ 
mal  conductivity.  The  Sutherland  Model  reads 


36.-  = 


(88) 


T  +  G*^  ’ 

where  tire  constants  and  G^  are  empirically  detramirred.  The  Gupta  Model  is  given  by 

%.  a  TlA*jOn'ri*+B3w(lnTP+C,jjlnT+D3„]  ^  (89) 


wlrere  A^,  C^j^,  and  are  the  curvefit  coefficients.  Again,  the  Sutherland  Model 

is  used  for  temperatures  below  1000  K. 


c.  MIXTURE  RULES 


Once  the  species  values  far  the  transport  properties  have  been  found,  then  an  appropriate 
mixing  rule  is  applied  to  give  the  total  mixture  property.  Two  mixture  rules  are  implemented  and 
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can  be  used  with  either  the  Sutherland  Model  or  Gupta  Model  for  the  transport  properties.  The 
first  mixing  rule  is  the  “standard”  Wilke’s  Rule  as  follows 

N 

(90) 

i=l 

where  the  weighting  function  is  given  by 


w, - ^ 


»  N 


j=l 


with  the  coefficient  <t>y  given  by  the  following  expression 


(91) 


In  the  above,  ^  ddnof  ^a  ansport  property  and  stands  for  either  viscosity  ot  thermal  conductiv¬ 
ity.  Also,  the  species  concentration  Xj  =  has  been  introduced.  This  mixing  rule  will  be 
referred  to  as  Wilke’s  Rule. 


Gupta  et  al.  [Reference  17],  prq)ose  an  improved,  but  more  complicated  mixing  rule  based 
on  the  collision  cross-sections.  The  mixture  rule  will  have  the  same  form  as  in  Equation  (90)  but 
the  weighting  functitm  is  now  given  by 


W,. 


N  yi  + 


j-1 


(93) 


where 


(94) 

Here,  the  collision  cross-sections  have  been  introduced  and  Aq^,  Bq^,  Cq^  and  D  are 
the  curvefit  coefficioits.  This  will  be  referred  to  as  Gupta’s  Rule. 
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SECTION  IV 


NUMERICAL  FORMULATION:  BLACK  BOX 


In  the  following,  details  of  the  numerical  algorithms  used  to  determine  the  equilibrium  com¬ 
position  are  given.  In  addition,  the  methods  developed  in  Section  n  to  obtain  the  thermodynamic 
properties  of  a  gas  mixture  in  local  chemical  equilibrium  are  further  illustrated.  The  possibility 
of  multiple  solutions  to  the  chemical  equilibrium  problem  and  the  correction  strategies  employed 
to  prevent  non-physical  results  are  also  discussed.  Finally,  methods  used  to  enhance  the  robust¬ 
ness  and  efficiency  of  the  “Black  Box”  are  detailed. 

1 .  “BLACK  BOX”  SOLVER 

The  first  step  in  the  numerical  solution  of  the  systems  of  nonlinear  govoning  equations  de¬ 
veloped  in  Section  n  involves  the  linearizadop.  of  the  equations  using  the  Newton-Raphson  tech¬ 
nique.  The  hnear  system  will  have  the  form 

Ax  =  b  ,  (95) 

where  A  is  the  (H  +  1)  X  (H  +  1)  Jacobian  matrix.  For  the  Mass  Constraint  Technique  H 
equals  N,  and  for  the  Degree  of  Advancement  Technique  H  will  be  equal  to  N  —  NE.  The  se¬ 
cond  step  comprises  the  direct  inversion  of  the  linear  system  by  means  of  a  LU  decomposition, 
where  partial  pivoting  is  employed.  Relative  and  absolute  reduction  of  the  residual  vector  is 
checked  and  iterations  are  performed  using  double-byte  arithmetic  until  either  has  reached  a  pre¬ 
scribed  tolerance.  For  use  with  the  flow  solver,  the  matrix  of  computational  cells  are  converted 
to  a  vector,  so  that  the  Black  Box  solves  one  system  of  equations  per  cell  in  the  flow  field,  and 
vectorizes  easily  over  the  number  of  computational  cells. 

a.  MASS  CONSTRAINT  TECHNIQUE 

The  specific  components  of  the  linearized  system  of  equations.  Equation  (95),  for  the  Mass 
Constraint  Technique,  written  for  constant  density  and  internal  energy,  will  have  the  following 
form 
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^^N-NE 

aQi 

^^N-NE+1 

aoi 


^^N-NE  ^^N-NE 

apN  aT 

^^N-NE+1  ^^N-NE+l 

^On 


=  [AQi,  ...,  Aqn.  at]  ,  (97) 

bT  =  -  [wi,  ....Wn-ne,Wn-ne+i . WN,d]  .  (98) 

In  the  above,  w  denotes  the  laws  of  mass  action  aial  the  mass  constraints,  and  d  denotes  the  enm*- 
gy  equation.  These  have  been  already  defined  by  Equations  (30-32).  In  addition,  the  partial  de¬ 
rivatives  of  the  laws  of  mass  action  will  read 


,  yu  nr  N  /  ^  N  /  ^ 

v-nfe  V. 

J  Jr=l  m=l'  '  m=l'  ' 


dT 


i  =  1 . N-  ME  , 


(100) 


where  in  Equation  (99) 

y'm^  (Vmj)  if  m  j  , 

vJnj  (viir)  =  '  vinj  ~  1  (vj^  -  1)  if  m  =  j  and  vjn,  (vi;,,)  0  , 

0  if  via,  (v^,)  =  0  . 

» 

The  partial  derivatives  of  the  mass  constraints  are 

^^N-NE+i  _  ^ 

=  0  ,  i=l . NE. 


(101) 


(102) 
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Finally,  the  partial  derivatives  of  the  equation  of  state  are  given  by 


I  Cv/t)dt  +  hf.  , 


(103) 

(104) 


b.  DEGREE  OF  ADVANCEMENT  TECHNIQUE 


The  specific  components  of  the  hneaiized  system  of  equations.  Equation  (95),  for  the  De¬ 
gree  of  Advancement  Technique,  written  for  constant  density  and  internal  energy,  will  have  the 
following  form 


dWj^ 

dWi 

dWi 

^^N-NE 

~dT 

1 

t 

\ 

(105) 

^^N-NE 

BT 

dd 

dd 

Bd 

all 

^^N-NE 

BT 

II 

• 

AT]  , 

(106) 

bT  =  _  [^1,  .. 

(107) 

The  advantage  of  the  Degree  of  Advancement  Technique  should  be  readily  apparent  in  the 
above.  The  system  has  been  reduced  from  N  +  1  to  N  -  NE  +  1  equations.  However,  while  the 
Jacobian  is  smaller  component-wise,  it  has  been  rendered  computationally  more  complex.  Using 
the  diffmential  form  of  the  transformation  relation.  Equation  (37),  the  partial  derivatives  with 
respect  to  the  degrees  of  advancement  can  be  written  in  terms  of  the  previous  entries  as  follows 


ni»l 


(108) 


|i  -  e  II  -  vinj)^  .  (109) 

1  in*  1 

where  the  partials  with  respect  to  temperature  remain  the  same  and  the  mass  constraint  equations 
have  been  eliminated. 
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2.  ALTERNATE  THERMODYNAMIC  SYSTEMS 


The  preceding  subsection  dealt  with  the  general  situation  of  the  complete  thermodynamic 
system  being  defined  by  density  and  internal  energy,  CTS  =  f(e,e).  As  was  discussed  in  Sec¬ 
tion  n,  situations  arise  where  the  complete  thermodynamic  system  is  known  as  a  function  of 
pressure  and  density,  CTS  =  f(p,  q),  or  as  a  function  of  pressure  and  temperature, 

CTS  =  f(p,  T).  In  the  following  discussion  on  these  alternate  thermodynamic  systems,  only  the 
Mass  Constraint  Technique  will  be  considered. 

For  the  situation  where  CTS  =  f(p,Q),  the  linearized  system  for  the  Mass  Constraint  Tech¬ 
nique,  given  by  Equations  (96-98),  remains  the  same,  with  the  exception  that  closure  is  now 
provided  by  the  thermal  equation  of  state 


N 

d(ej,T)  =  Y,  OjRjT  -  (p)o  =  0 ,  (110) 

j=i 

where  (p)o  is  the  initial  pressure.  The  partial  derivatives  of  the  equation  of  state,  previously  giv¬ 
en  by  Equations  (103)  and  (104),  are  now  written  as 


RjT 


(111) 

(112) 


The  situation  where  pressure  and  temperatiire  are  known,  CTS  =  f(p,  T),  requires  more 
modification.  As  stated  previously,  the  dependent  variables  are  now  species  densities  and  total 
density,  and  the  thermal  equation  of  state  has  replaced  the  caloric  equation  of  state  in  the  govern¬ 
ing  equations.  Although  total  density  is  nothing  more  than  the  sum  of  the  species  densities  and 
would  indicate  that  one  equation  could  be  dropped  from  the  system,  reducing  the  size  of  the 
problem  is  not  practical.  Li  fact,  the  state  relationship  is  required  to  deltine  the  thermodynamic 
state  and  while  it  could  be  substituted  into  the  rate  equations,  this  would  increase  the  complexity 
of  the  Jacobians,  making  them  more  expensive  to  compute.  The  linearized  system  for  constant 
pressure  and  temperature  will  be 
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(113) 
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(114) 


and  the  residual  vector  remains  the  same.  The  partial  derivatives  of  the  laws  of  mass  action  with 
respect  to  the  species  densities  remain  the  same  and  are  given  by  Equation  (99),  whereas  the  par¬ 
tial  derivatives  of  the  laws  of  mass  action  with  respect  to  total  density  are  zero.  The  partial  de¬ 
rivatives  of  the  mass  constraints  are  given  as 


^^N-NE  +  i  _  ^ 


(115) 


^^N-NE+i  _ 
dQ 


i  =  1, 


NE. 


(116) 


Recognizing  that  pressure  and  temperature  are  constant,  the  thermal  equation  of  state  can  be 
written  in  the  form 


N 

<i(Qj,e)  =  “  (p)o  =  0 »  (117) 

where  the  terms  denoted  by  (  )o  rqnesent  initial  conditions.  The  partial  derivatives  of  the  ther¬ 
mal  equation  of  state  will  then  be 


(118) 


(119) 
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3.  THERMODYNAMIC  PROPERTIES 


Once  the  composition  of  a  mixture  of  thermally  perfect  gases  has  been  evaluated,  the  ther¬ 
modynamic  properties  of  the  mixture  can  be  determined.  These  are  primarily  functions  of  the 
temperature  and  density  derivatives  of  the  species  densities.  As  was  discussed  in  Section  II,  the 
determination  of  these  derivatives  is  a  straightforward  process  involving  the  solution  of  two  lin¬ 
ear  systems  of  algebraic  equations,  provided  that  temperature  and  density  derivatives  of  the  gov¬ 
erning  equations  are  readily  available.  This  is  indeed  the  case  when  the  Newton-Raphson  algo¬ 
rithm  is  utilized. 

Taking  the  Mass  Constraint  Technique  as  an  example,  it  should  be  readily  apparent  that  the 
first  N  rows  and  columns  of  the  converged  Jacobian  matrix  corresponds  to  the  matrix  A  required 
to  determine  the  partial  derivatives.  Also,  the  last  column  of  the  Jacobian  matrix  is  the  RHS,  giv¬ 
en  by  C  in  Equation  (47),  for  the  temperature  derivative  computation.  The  RHS  for  the  density 
dervatives,  given  by  B  in  Equation  (46),  is  also  easily  derived  and  will  be 

^=0,  i=l . N-NE,  (120) 

and 


'^’^N-NE+i 

dQ 


i  =  1 . NE.  (121) 


Similar  considerations  can  be  ^plied  to  tlte  Degree  of  Advancement  Technique,  where  the 
matrix  A  is  now  given  by  the  first  N  -  NE  rows  and  columns  of  the  Jacobian  matrix  at  conver¬ 
gence.  The  RHS  terms  will  be  the  same  as  those  for  the  Mass  Constraint  Technique.  Therefore, 
at  convergence,  the  temperature  and  density  derivatives  are  easily  evaluated.  At  this  point,  the 
reason  behind  the  choice  of  density  and  temperature  as  the  two  fundamental  state  variables  to 
describe  the  properties  of  a  mixture  in  local  chemical  equilibrium  becomes  apparent.  If  one  were 
to  choose  density  and  internal  energy,  for  example,  evaluation  of  the  partial  derivatives  with  re¬ 
spect  to  e  would  be  computationally  more  expensive  than  the  method  presented  above,  due  to  the 
fact  that  the  laws  of  mass  action  are  not  explicitly  defined  in  terms  of  internal  energy. 


4.  MULTIPLE  SOLUTIONS 


llie  possibility  of  multiple  solutions  is  an  intoesting  yet  troublesome  problem  that  affects 
chemical  equilibrium  calculations.  The  nonlinear  nature  of  the  governing  equations  potentially 
allows  several  different  mathematical  solutions  for  the  same  set  of  equations. 

An  example  of  this  occurrence  can  be  given  for  the  5-Species  Air  Model  using  the  Vibra¬ 
tional  Model  and  Curvefit  Kc.  The  initial  guess  for  the  temperature  is  5000  K  and  the  initial 
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guesses  for  the  mass  fractions  are  the  exact  values,  where  density  q  =  8 . 676  x  10  ^kg/m^ 
and  internal  energy  e  =  6 . 271  x  lO^kJ/kg  have  been  specified.  Table  1  gives  the  temperature 
and  composition  for  two  solutions  obtained  for  the  problem  defined  above.  Although  both  math¬ 
ematically  satisfy  the  laws  of  mass  action,  mass  constraints  and  energy  equation,  the  first  solu¬ 
tion  is  physically  not  possible  since  mass  firactions  cannot  have  negative  values.  Physically,  there 
can  be  only  one  solution  with  nonnegative  mass  fi'actions,  since  the  specification  of  two  state 
variables  describes  a  unique  thermodynamic  state.  The  above  example  illustrates  the  possibility 
of  negative  densities,  but  it  is  also  possible  to  obtain  solutions  with  negative  temperatures. 

In  light  of  the  problem  mentioned  above,  three  techniques  were  develc^jed  to  guide  conver¬ 
gence  to  the  one  physical  solution.  The  first  two  methods  employ  a  scaling  factor  to  reduce  the 
magnitude  of  the  species  densities  updates,  AQ^  =  Qj*'*'  ^  -  Qp  where  the  superscript  “  rqsre- 
sents  the  iteration  level.  The  first  method  will  be  called  Catastrophic  Limiter  and  it  intervenes 
only  when  a  species  density  is  going  to  become  negative  from  the  iteration  update.  The  second 
method  limits  the  magnitude  of  the  relative  change  in  species  densities  at  every  iteration,  and  is 
called  Relative  Limiter.  If  the  same  scaling  factor  is  used  on  all  the  species  densities  updates, 
then  it  can  be  shown  that  the  mass  constraints  will  not  be  violated.  This  is  important  when  using 
these  methods  with  the  Degree  of  Advancement  Technique,  which  does  not  explicitly  enforce  the 
mass  constraints.  The  scaling  factor  a  will  have  the  form 


f 

1  P  e?  ' 
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> 

where  P  is  the  maximum  allowed  relative  change  in  species  densities 


(122) 


^^new 


i  =  1,...,N 


(123) 


The  only  difference  between  the  two  methods  is  that  Relative  Limita*  applies  the  correction  at 
each  iteration,  while  Catastrt^hic  Limiter  applies  the  correction  only  when  a  species  density 
would  become  negative  as  a  result  of  the  update.  In  order  to  prevent  the  convergence  process 
firom  slowing  down  dramatically,  the  scaling  factor  is  limited  to  a  minimum  value,  o^.  The 
new  value  for  the  species  densities  updates  will  then  be  given  by 


=  o  Aq?“  ,  i  =  1,...,N .  (124) 

Neith^  of  the  above  methods  can  prevent  negative  species  densities  from  occurring  in  die  tran¬ 
sient  when  extenuating  circumstances  occur  (i.e.  q"  =  0  and  AQj  <  0).  However,  the  negative 
species  densities  are  omitted  from  the  scaling  factor  evaluation,  see  Equation  (122),  since  en- 
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forcing  a  scaling  factor  on  the  negative  species  densities  would  prevent  those  densities  from  be¬ 
coming  positive  again. 

While  the  use  of  a  common  scaling  factor  ensures  satisfaction  of  the  mass  constraints,  the 
linearized  equation  of  state  would  not  be  satisfied  by  the  modified  species  densities  updates,  un¬ 
less  the  temperature  update  is  also  suitably  modified.  This  can  be  accomplished  by  resolving  the 
last  row  of  the  governing  system  for  the  modified  temperature  update,  using  the  modified  spe¬ 
cies  densities  updates. 

The  third  correction  method  is  based  on  the  absolute  Newton’s  method  described  by 
Meintjes  and  Morgan  [Reference  4],  and  will  be  called  Absolute  Newton  Limiter.  Negative  spe¬ 
cies  are  avoided  by  taking  the  new  values  of  tl%  species  densities,  to  be  the  absolute  value  of  the 
sum  of  the  old  values  and  their  updates,  q?'*'  ^  =  10"  +  AQj  | .  This  method  does  not  use  a 
common  scaling  factor  and  hence  violation  of  the  mass  constraints  is  possible  during  the  itera¬ 
tive  process.  Due  to  this  violation,  the  Absolute  Newton  Limiter  is  unsmtable  for  use  with  the 
Degree  of  Advancement  Technique,  as  will  be  shown  in  the  results. 

Another  occurrence  associated  with  the  Newton-R{q>hson  method  is  the  over-shooting  of 
the  temperature,  which  can  lead  to  negative  temperatures.  This  is  especially  possible  for  ill- 
posed  problems  involving  low  temperature  compositions  where  the  initial  guess  for  temperature 
is  given  as  a  large  value.  When  written  in  terms  of  temperature  instead  of  species  density,  the 
Relative  Limiter  and  Absolute  Newton  Limiter  can  be  used  to  prevent  these  over-shootings. 

5.  ROBUSTNESS  AND  EFnCffiNCY 

Several  techniques  were  employed  to  enhance  the  robustness  and  efficiency  of  the  Black 
Box.  These  techniques  include  improving  the  initial  estimates  of  the  dependent  variables,  reduc¬ 
ing  the  size  of  the  vector  of  computational  cells  solved  by  the  Black  Box,  and  hneezing  the  chem¬ 
istry. 


a.  INITIAL  GUESS 


One  of  the  key  factors  in  the  robustness  and  efficiency  of  a  numerical  scheme  employing  the 
Newton-Raphson  method  is  the  availability  of  **good”  initial  guesses.  Two  methods  are 
employed  in  the  Black  Box  to  provide  these  “good”  initial  guesses  for  the  species  densities  and 
temperature.  The  first  one  consists  of  assigning  the  values  obtained  at  the  previous  time  stq>  to 
tl»  new  values.  Thus  given  new  values  for  density  and  total  energy  per  unit  volume  from  the 
flow  solver,  the  initial  guesses  for  the  species  densities  and  temperature  will  be 


jvi+1  _ 
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where  the  superscript  "  again  denotes  the  time  step. 

The  second  method  uses  the  mass  fraction  derivatives  to  provide  an  improved  initial  guess. 
The  updated  values  will  be  given  by 
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6/e  \ 

'<?  . 

=  T"  +  (f  )^Ae  +  (f  )"Ae  .  (127) 

where  the  partial  derivatives  are 


and  Aq  =  q"'*'  ^  -  q"  and  Ae  =  e"'*'  ^  —  e"  are  provided  by  the  flow  solver.  The  partial  deriva¬ 
tives  of  the  species  mass  fractions  with  respect  to  density  and  temperature  are  obtained  from  the 
approach  already  described. 

b.  VECTOR  REDUCTION 

A  very  effective  method  used  to  increase  the  efficiency  of  the  Black  Box  is  vector  reduction, 
whereby  the  vector  of  computational  cells  is  reduced  in  size  as  cells  converge.  Considering 
three-dimensioiud  flow  simulations,  initially  all  the  computational  cells  are  stored  in  an  active 
vector  with  a  vector  length  equal  to  the  product  of  the  three  dimensions  of  the  computational 
matrix  imder  consideration.  The  residuals  for  each  of  the  governing  systems  is  checked  at  each 
iteration  for  convergence.  As  a  cell  converges  it  is  removed  from  the  active  vector  and  the  vectm* 
length  is  correspondingly  reduced.  Thus  the  expense  associated  with  the  cmnputation  of  the  Ja- 
cobians  for  converged  cells  is  not  incurred  aixl  iterations  are  performed  on  the  imconverged  cells 
only.  The  efficiency  of  vector  reduction  does  not  come  without  cost.  In  order  to  implement  the 
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vector  reduction,  a  memory  intensive  data  structure  is  used  to  keep  track  of  the  cells  in  the  vec¬ 
tor  and  their  corresponding  positions  in  the  computational  matrix.  However,  it  will  be  seen  in  the 
results  section  that  the  benefits  definitely  outweigh  the  costs. 

c.  FREEZING  THE  CHEMISTRY 


As  should  be  {^parent,  use  of  the  Black  Box  to  provide  the  composition  and  thermodynamic 
properties  of  a  gas  mixture  for  use  in  a  flow  solver  requires  some  computational  expense.  This 
expense  can  be  reduced  by  freezing  the  chemistry  for  a  desired  number  of  itmtions  of  the  flow 
solver,  where  the  chemical  composition  would  be  evaluated  at  every  k^**  iteration  instead  of  at 
every  iteration  of  the  flow  solver.  At  the  iterations  where  the  chemistry  is  frozen,  the  mass  frac¬ 
tions  are  considered  to  be  constants.  Thus  the  governing  system  of  equations  in  the  Black  Box 
would  reduce  to  an  equation  of  state  only  and  the  unknowns  would  reduce  to  just  temperature  or 
total  density,  depending  on  the  known  state  variables.  Iterations  still  need  to  be  performed  for  die 
general  case,  CTS  =  f(Q,  e),  since  the  caloric  equation  of  state  is  a  nonlinear  function  of  temper¬ 
ature  and  cannot  be  solved  diiecdy.  However,  solutions  for  the  cases  CTS  =  f(p,  q)  and 
CTS  =  f(p,  T)  are  easily  obtained,  where  the  thermal  equation  of  state  can  be  solved  direcdy.  In 
the  first  case,  temperature  is  readily  given  as 


and  for  the  second  case,  density  is  easily  obtained  by 


(132) 


(133) 


where  the  light-haiKl  sides  of  both  equaticxis  are  known  constants.  In  practice,  the  number  of 
iterations  at  frozen  chemistry  will  be  small  at  the  beginning  of  the  flow  computation,  when  the 
thermodynamic  state  and  hence  chemical  composition  is  changing  rapidly.  Towards  cmiver- 
gence,  this  number  can  be  significantly  increased,  and  the  overall  computational  time  reduced. 
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Temp.  (K) 

N2  Mass  Fraction 
O2  Mass  Fraction 
NO  Mass  Fraction 


N  Mass  Fraction 


O  Mass  Fraction 


Negative 

3932.45 

0.803368 

0.079529 

-0.076377 

-0.000622 

0.194102 


Exact 

4000.00 

0.740236 

0.042519 

0.055901 

0.000761 

0.160582 


Table  1  Example  of  multiple  solution. 
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SECTION  V 


BLACK  BOX  RESULTS 

Numerical  results  were  obtained  for  each  of  the  chemistry  models  detailed  in  Section  Ed  us¬ 
ing  different  densities  over  varying  ranges  of  temperature.  Unless  otherwise  stated,  all  calcula¬ 
tions  were  performed  using  the  Vibrational  Model  and  Consistent  Kc. 

1.  AIR  MODELS 

In  the  following,  equilibrium  compositions  for  air  were  computed  at  densities  q  =  1.293, 
1.293  10“^  and  1.293  10“^  kg/m^,  where  the  first  value  corresponds  to  standard  conditions. 
The  compositions  computed  using  S-Species  Air  Model,  plotted  as  mole  percentage  versus  tem¬ 
perature,  are  shown  in  Figure  1,  Figure  2  and  Figure  3  fm*  the  three  respective  densities.  Of  im¬ 
portance  is  the  result  that  as  density  decreases,  so  does  the  temperature  at  which  the  mixture  has 
totally  dissociated  into  monatomic  components.  For  the  density  q  =  1.293  kg/m^,  the  mixture 
still  contains  diatomic  species  up  to  16,000  K,  as  seen  in  Figure  1.  The  equilibrium  composition 
for  density  e  =  1.293  10  kg/m^,  as  seen  in  Figure  2,  becomes  almost  completely  monatom¬ 
ic  at  around  14,000  K.  Also,  the  mole  percentage  of  NO  peaks  at  around  3,000  K,  which  is  con¬ 
sistent  with  the  findings  of  other  researchers  [Reference  23].  In  Figure  3,  where  the  composition 
for  density  Q  =  1.293  10  kg/m^  is  plotted,  it  can  be  sern  that  the  temperature  at  which  the 
mixture  is  virtually  dissociated  has  dropped  to  around  9,000  K. 

The  variation  of  traiperature  versus  internal  energy,  computed  using  5-Species  Air  Model, 
is  shown  in  Figure  4  fen*  the  tinee  density  values  given.  The  effect  of  chemical  reactions  are 
readily  iq^parent  in  this  figure.  For  an  ideal  gas,  temperature  would  be  a  linear  fimetion  of  inter¬ 
nal  energy  and  independent  of  density. 

The  variation  of  isentropic  itxiex  versus  temperature  is  given  in  Figure  5  for  the  same  three 
densities  and  same  air  model.  At  low  temperatures,  the  index  starts  at  its  diatomic  limit  of  7/5.  It 
then  follows  a  series  of  peaks  and  valleys  until  it  finally  reaches  its  monatomic  limit  of  5/3. 
Comparing  each  dmsity  curve  with  its  corresponding  composition  plot  given  in  Figure  1  through 
Figure  3,  it  can  be  seen  that  the  first  valley  correlates  well  with  the  point  of  massive  O2  dissoci- 
atioiL  The  subsequent  peak  ccnrelates  with  a  peak  in  O  production.  The  second  valley  correlates 
with  the  point  of  massive  N2  dissociation.  It  can  be  inferred  firom  Figure  5  that  it  may  be  hazard¬ 
ous  to  assume  a  constant  “gamma”  in  the  derivatiem  of  flux-^lit  algoritiims.  The  variation  of 
speed  of  sound  versus  temperature  is  shown  in  Figure  6.  Again,  deviation  fixnn  itteal  behavior  is 
noticeable,  as  this  plot  would  reduce  to  a  single  parabolic  shape  for  an  ideal  gas. 
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Equilibrium  compositions  for  density  Q  =  1.293  10~^  kg/m^  are  shown  in 
Figure  7  through  Figure  10  for  the  7, 9,  11  and  13-Species  Air  Model  respectively.  The  7-Spe¬ 
cies  Air  Model  solution  given  in  Figure  7  is  nearly  identical  to  Figure  2  except  for  the  two  addi¬ 
tional  species  NO  and  e  “ .  The  mole  percent  for  these  two  species  is  identical  due  to  preserva¬ 
tion  of  electrical  chaise.  However,  it  should  be  noted  that  the  physical  behavior  of  the  7-Species 
mixture  is  inaccurate,  since  at  the  high  end  of  the  temperature  range  the  mixture  should  be  en¬ 
tirely  monatomic.  This  is  because  the  7-Species  Air  Model  has  no  mechanism  for  the  production 
of  atomic  ions,  specifically  N  and  O'*” .  A  more  noticeable  change,  other  than  the  addition  of 
more  species,  occurs  in  the  more  complex  models  which  incorporate  reactions  describing  the 
production  of  atomic  ions.  At  the  higher  temperatures,  the  mixture  is  composed  entirely  of  mo¬ 
natomic  species,  as  can  be  seen  in  Figure  8  through  Figure  10  .  Furthennore,  these  mixtures  are 
becoming  increasingly  more  ionized. 

Plots  of  mole  percentage  versus  temperature  for  the  17-Species  Air  Model  and  the  three 
densities  given  previously  are  depicted  in  Figure  11 ,  Figure  12  and  Figure  13.  Again  the  depen¬ 
dence  of  the  level  of  dissociation  on  density  is  readily  apparent.  These  plots  are  in  excellent 
agreement  with  similar  results  obtained  by  Hilsenrath,  Klein  and  Wooley  [Reference  24],  where 
the  equilibrium  computations  were  performed  using  27  species. 

The  variation  of  temperature  versus  internal  energy  is  given  in  Figure  14  for  the  three  densi¬ 
ties  and  17-Species  Air  Model.  The  curves  are  nearly  identical  to  those  computed  using  5-Spe¬ 
cies  Air  Model,  shown  in  Figure  4,  up  to  an  internal  energy  e  =  40  MJ/kg.  At  this  value  for  inter¬ 
nal  energy  the  curves  for  the  5-Species  Air  Model  reach  their  limiting  linear  functional  form, 
since  the  composition  is  fairly  constant  past  that  point.  However,  for  the  17-Species  Air  Model, 
this  limit  is  not  reached  until  around  200  MJ/kg.  The  reason  for  this  difference  is  that  the  more 
complex  air  model  incorporates  charged  species,  where  internal  energy  is  used  in  the  ionization 
reactions. 

Plots  of  the  isentropic  index  versus  temperature  for  the  same  three  values  of  density  and 
17-Species  Air  Model  are  given  in  Figure  15.  As  in  the  results  for  5-Species  Air  Model,  the 
isentropic  index  starts  at  its  diatomic  limit  and  follows  a  series  of  peaks  and  valleys  until  it 
reaches  its  monatomic  limit.  Again,  when  compared  to  the  corresponding  plots  of  composition 
given  by  Figure  11  through  Figure  13,  the  first  two  valleys  and  first  peak  correlate  in  a  similar 
manner  as  the  results  for  the  5-Species  Air  Model.  Additionally,  the  second  peak  correlates  vrith 
the  peak  in  N  production  and  the  subsequent  valley  correlates  nicely  with  the  beginning  of  mas¬ 
sive  ionization.  The  variation  of  speed  of  sound  versus  temperature  for  the  17-Species  Air  Mod¬ 
el  is  shown  in  Figure  16. 
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Comparisons  of  the  air  models  are  made  in  Figure  17  and  Figuie  18  for  density 
Q  =  1.293  10“^  kg/m^.  The  variation  of  temperature  versus  internal  energy  for  all  the  air 
models  is  given  in  Figure  17.  Here  the  effects  of  ionization  incorporated  in  the  more  complex 
models,  9,  11,  13  and  17-Species  Air  Model,  are  readily  apparent.  The  variation  of  isentropic 
index  versus  temperature  is  shown  in  Figure  18.  The  inaccuracy  of  the  7-Species  Air  Model  is 
very  apparent  in  this  plot,  where  at  high  temperatures  the  isentropic  index  tends  towards  the  di¬ 
atomic  limit  and  then  drops  off.  Moreover,  it  should  be  pointed  out  that  for  the  simpler  air  mod¬ 
els  the  diatomic  limit  is  reached  at  a  much  lower  temperature  than  for  the  more  complex  models. 

It  is  important  to  point  out  that  the  results  from  all  the  chemistry  models  are  the  same  up  to 
a  temperature  corresponding  to  the  onset  of  ionization,  which  for  the  density 
Q  =  1.293  10“^  kg/m^  results  shown  is  about  8,000  K.  At  higher  temperatures,  where  ioniza¬ 
tion  effects  are  important,  the  9, 11, 13  and  17-Species  Air  Models  are  in  perfect  agreement. 
Thus  for  high-speed,  high-temperature  air  flow  simulations,  there  should  be  no  loss  in  accuracy 
through  the  use  of  the  simpler  9-Species  Air  Model  instead  of  the  computationally  more  expen¬ 
sive  17-Species  Air  Model,  Moreover,  if  the  maximum  temperature  of  the  flow  field  is  less  than 
the  temperature  at  which  ionization  occurs,  then  further  computational  savings  can  be  realized 
through  the  use  of  the  simple  5-Species  Air  Model  without  any  appreciable  loss  in  accuracy. 

2.  COMBUSTION  MODELS 

The  equilibrium  composition  obtained  for  Hydrogen- Air  Combustion  Model  is  given  in 
Figure  19  for  a  density  q  =  1.324  10“^  kg/m^  The  equivalence  ratio  ^  =  0.29841  corre¬ 
sponds  to  a  lean  mixture.  The  accuracy  of  this  model  is  questionable,  since  at  the  higher  temper¬ 
atures  oxygen  should  be  dissociated. 

The  variation  of  temperature  versus  internal  energy  for  Hydrogen-Air  Combustion  Model  is 
shown  in  Figure  20  for  densities  of  Q  =  1.324  10“  ^  1.324  10“^  and  1.324  10“^  kg/m^.  The 
variation  of  speed  of  sound  versus  temperature  is  given  in  Figure  21  for  the  same  three  densities. 
It  can  be  seen  that  the  thermodynamic  properties  of  the  mixture  are  not  strongly  affected  by  den¬ 
sity  for  this  combustion  model.  This  is  probably  due  to  the  high  concentration  of  inert  nitrogen 
in  the  mixture. 

Equilibrium  compositions  computed  for  Hydrogen-Oxygen  Combustion  Model  are  shown 

in  Figure  22  and  Figure  23  for  densities  Q  =  2.8306  and  2.8306  10“^  kg/m^,  and  a  mixture 
ratio  of  6: 1 .  This  is  the  same  composition  used  in  the  space  shuttle  main  engine  (SSME)  nozzle 
studies  of  Wang  and  Chen  [Reference  13].  At  low  temperatures,  water  and  excess  diatomic  hy¬ 
drogen  are  the  prevalent  species,  while  at  higher  temperatures,  where  the  dissociation  reactions 
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are  predominant,  atomic  hydrogen  and  oxygen  are  noticeable.  This  figure  illustrates  the  depen¬ 
dence  on  temperature  as  to  whether  the  mixture  is  fully  or  partially  combusted. 

The  variation  of  temperature  versus  internal  energy  for  the  Hydrogen-Oxygen  Combustion 
Model  is  shown  in  Figure  24  for  the  two  densities  given  above.  A  much  stronger  dependence  on 
density  is  registered  for  this  model  as  opposed  to  the  previous  model.  The  variation  of  isentropic 
index  versus  temperature  is  given  in  Figure  25  for  the  same  two  densities  and  same  combustion 
model.  When  compared  to  the  composition  plots  of  Figure  22  and  Figure  23,  the  minimas  of  the 
isentropic  index  curves  correlate  with  the  points  where  oxygen  is  beginning  to  massively  dissoci¬ 
ate.  The  composition  of  both  mixtures  are  fairly  constant  up  to  a  temperature  around  3,000  K. 
This  is  reflected  in  the  plot  of  speed  of  sound  versus  temperature,  given  in  Figure  26,  where  the 
curves  for  the  two  densities  begin  to  deviate  at  this  temperature. 

The  equilibrium  composition  obtained  using  Hydrocarbon  Combustion  Model  for  density 
Q  =  3.295  10”^  kg/m^  is  shown  in  Figure  27,  where  the  equivalence  ratio  <I>  =  1.2  corre¬ 
sponds  to  a  fuel-rich  mixture.  This  is  the  same  composition  studied  by  Meinqes  and  Morgan 
[Reference  4].  Again,  the  dependence  of  the  combustion  process  from  the  temperature  is  readily 
t^parent. 

The  variation  of  temperature  with  internal  energy  for  the  Hydrocarbon  Combustion  Model 
is  shown  in  Figure  28,  where  three  curves  are  plotted  for  densities  q  =  3.295,  3.295  10“^  and 
3.295  10“^  kg/m^.  Similar  to  the  comparisons  made  earlier,  the  valley  in  the  density 
Q  =  3.295  10"^  kg/m^  curve,  in  the  plot  of  isentropic  index  versus  temperature  shown  in  Fig¬ 
ure  29,  is  seen  to  correlate  well  to  the  point  of  massive  dissociation  of  O2  seen  in  Figure  27.  The 
variation  of  speed  of  sound  is  given  in  Figure  30  for  the  same  three  densities. 

3.  PLASMA  MODEL 

Equilibrium  compositions  obtained  from  the  Argon  Plasma  Model  are  plotted  in  Figure  31 
for  the  densities  Q  =  1.303  10“^  1.303  10“^  and  1.303  10“^  kg/m^  proceeding  from  left 
to  right.  Plots  of  the  variation  of  temperature  versus  internal  energy,  isentropic  index  versus  tem¬ 
perature  and  speed  of  sound  versus  temperature  are  given  in  Figure  32,  Figure  33  and  Figine  34 
respectively,  for  the  same  three  values  of  density.  A  strong  dependence  on  density  is  registered  in 
all  the  plots. 

4.  TRANSPORT  PROPERTIES 

Comparisons  of  the  transport  property  evaluation  methods  were  made  using  the  1 1 -Species 
Air  Model  under  a  constant  pressure  of  one  atmosphere.  Either  of  the  two  models  used  to  com- 
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pute  the  species  transport  property,  Sutherland  Model  or  Gupta  Model,  can  be  combined  with 
either  of  the  mixture  rules,  Wilke  Rule  or  Gupta  Rule,  to  give  the  value  for  the  mixture.  The 
variation  of  the  mixture  viscosity  versus  temperature  is  given  in  Figure  35  for  each  of  four  per¬ 
mutations  possible,  while  the  variation  of  mixture  thermal  conductivity  is  given  in  Figure  36.  In 
both  plots  the  values  have  been  normalized  using  the  appropriate  Sutherland’s  Law  for  air  as  a 
homogeneous  mixture.  There  is  relatively  little  difference  between  the  models  at  low  tempera¬ 
tures.  However,  at  the  higher  temperatures,  the  Gupta-Gupta  technique,  which  incoiporates  col¬ 
lision  cross  sections,  could  prove  to  be  more  accurate. 

5.  ROBUSTNESS 

The  problem  of  multiple  solutions  to  the  equilibrium  equations  and  the  limiters  developed  to 
correct  this  problem  were  discussed  in  Section  IV.  Five  different  strategies  were  developed,  in¬ 
corporating  the  limiters  in  various  combinations  in  order  to  improve  the  robustness  of  the  Black 
Box.  These  correction  strategies  were  tested  using  Ideal  Dissociating  Oxygen  Model  as  weU  as 
the  5  and  11-Species  Air  Models.  Initial  guesses  for  species  densities  and  temperature  for  each 
cell  were  selected  with  a  certain  arbitrariness  and  all  calculations  were  performed  driving  the  re¬ 
sidual  to  machine  zero,  using  double  byte  arithmetic. 

In  the  following.  Strategy  1  will  refer  to  the  baseline  solution,  where  no  limiters  were  used. 
Strategy  2  uses  the  Catastrophic  Limiter  for  density  correction  and  the  Relative  Limiter  for  tem¬ 
perature  correction.  The  Relative  Limiter  was  used  for  both  density  and  temperature  correction 
in  Strategy  3.  Strategy  4  uses  Absolute  Newton  Limiter  for  both  density  and  temperature  correc¬ 
tion.  Strategy  5  is  similar,  except  temperature  correction  is  performed  using  Relative  Limiter. 
Preliminary  numerical  experiments  were  perfoimed  on  the  value  of  P  used  in  Equation  (122) 
and  the  value  P  =  1/3  was  found  to  be  the  optimal  choice. 

a.  OXYGEN  DISSOCIATION  RESULTS 

The  Ideal  Dissociating  Oxygen  Model  was  utilized  to  evaluate  the  equilibrium  composition 
at  three  diffoent  thermodynamic  states,  where  each  is  characterized  by  a  density 

e  =  3.0  10“^  kg/m^.  Case  1  corresponds  to  heavy  dissociation  and  intenuQ  energy 
e  =  8.14398  MJ/kg.  The  equilibrium  temperature  and  degree  of  dissociation  are 
T  =  2775.15  K  and  a  =  0.386874,  respectively.  Case  2  corresponds  to  light  dissociation  and 
internal  energy  e  =  3.0758  MJ/kg,  where  the  equilibrium  temperature  and  degree  of  dissoci¬ 
ation  are  T  =  2377.85  K  and  a  =  0.079069,  respectively.  Case  3  corresponds  to  almost  total 
dissociation  and  internal  energy  e  =  17.0  MJ/kg.  The  equilibrium  temperature  and  degree  of 
dissociation  are  T  =  3392.0  K  and  a  =  0.928500,  respectively. 
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The  results  of  the  robustness  study  for  the  Ideal  Dissociating  Oxygen  Model  are  presented  in 
Table  2  and  Table  3  for  the  Mass  Constraint  Technique  and  Degree  of  Ad'vapreme.  i:  Technique, 
respectively.  For  each  category,  the  percentage  of  results  that  are  correct,  blowup,  wrong  and  ex 
ceed  the  maximum  number  of  iterations  (50)  are  given,  and  are  denoted  in  the  table  by  %C,  %B, 
%W  and  %M,  respectively.  A  total  of  10,1(X3  computations  were  performcu  tor  each  case,  strate¬ 
gy  and  technique.  The  initial  guesses  for  temperature  and  degree  of  dissociation  were  varied  in 
the  range  100  K  <  T  ^  10,000  K  and  0  ^  a  5  1  in  equal  increments  of  100  K  for  the  tem¬ 
perature  and  0.01  for  the  degree  of  dissociation.  It  should  be  ^parent  that  Case  2  is  the  most  dif¬ 
ficult  case  of  the  three  since  it  has  the  lowest  percentage  of  solutions  that  are  correct  when  no 
coixections  are  made.  Strategy  2  is  completely  successful  at  guiding  convergence  to  the  correct 
solution.  Strategies  4  and  5  woric  perfectly  with  the  Mass  Constraint  Technique,  but  their  perfor¬ 
mance  with  the  Degree  of  Advancement  is  worse  than  no  correction  at  all.  This  was  expected, 
since  the  Degree  of  Advancement  does  not  enforce  the  mass  constraints  directly.  The  result  is 
solutions  where  density  is  not  conserved.  Strategy  3  is  not  quite  as  successful  as  Strategy  2.  The 
limit  on  the  relative  change  per  iteration  slows  convergence  down,  especially  for  particularly 
poor  initial  estimates  for  species  densities.  Hence,  Strategy  3  shows  an  increase  in  computations 
where  the  maximum  number  of  iterations  have  been  exceeded,  as  is  seen  in  the  table.  Also 
shown  in  Table  2  is  the  average  number  of  iterations  for  the  computations  converging  to  the  cor-  , 
rect  solution,  which  rarely  exceeds  12.  As  was  expected  Strategy  3  appears  to  retard  conver¬ 
gence.  Strategies  2  and  5  seem  to  be  the  most  efficient  methods. 

b.  AIR  DISSOCIATION  AND  IONIZATION  RESULTS 

In  a  similar  fashion,  robustness  studies  were  peifonned  using  the  5  and  1 1 -Species  Air 
Models.  Three  cases  were  used,  charactenzed  by  a  density  Q  =  1.293  10”^  kg/m^.  Case  1  cor¬ 
responds  to  massive  oxygen  dissociation  and  internal  energy  e  =  6.94  MJ/kg,  where  the  equi¬ 
librium  temperature  is  4000  K.  Case  2  corresponds  to  a  halfway  complete  nitrogen  dissociation 
and  internal  energy  e  =’  26.65  MJ/kg,  where  the  equilibrium  temperature  is  7(XX)  K.  Complete 
dissociation  of  nitrogen,  disiqipearance  of  NO  and  appearance  of  some  chaiged  species  charac¬ 
terize  Case  3,  where  internal  energy  e  =  36.621  MJ/kg,  and  the  equilibrium  temperature  is 
9000  K. 

Initial  guesses  for  temperature  were  varied  in  equal  increments  of  300  K  over  the  range 
300  K  s  T  s  30,(X)0  K.  Five  sets  of  initial  compositions  were  used:  the  fiorst  three  are  the 
equilibrium  compositions  corresponding  to  the  three  cases  above;  and  toe  final  two  compositions 
correspond  to  totally  undissociated  and  dissociated  nitrogen  and  oxygen  states,  respectively. 
These  initial  compositions  are  given  in  two  sets  denoted  Set  1  and  Set  2,  comprising  toe  first 
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three  compositions  and  all  five  compositions,  respectively.  Set  1  has  a  total  of  300  computations, 
while  Set  2  encompasses  all  500  computations.  In  light  of  the  Ideal  Dissociating  Oxygen  results, 
computations  were  performed  using  only  Strategies  2  and  5  along  with  the  baseline  Strategy  1. 

Results  for  5-Species  Air  Model  are  shown  in  Table  4  and  Table  5  for  the  Mass  Constraint 
Technique  and  Degree  of  Advancement  Technique,  respectively.  Notably  apparent  is  the  lack  of 
robustness  when  no  correction  strategies  are  implemented,  especially  for  Case  1.  Strategy  2  pro¬ 
vides  a  significant  improvement  over  the  baseline,  particularly  when  used  with  the  Degree  of 
Advancement  Technique.  Strategy  5  is  completely  successful  when  used  in  conjunction  with  the 
Mass  Constraint  Technique.  However,  as  was  already  seen  in  the  oxygen  dissociation  results,  it  is 
ill-suited  for  use  with  the  Degree  of  Advancement  Technique. 

Results  obtained  using  11-Species  Air  Model  are  given  in  Table  6  and  Table  7  for  the  Mass 
Constraint  Technique  and  Degree  of  Advancement  Technique,  respectively,  and  follow  the  same 
trend  as  the  previous  results  for  the  5-Species  Air  Model.  However,  the  overall  percentage  of 
correct  results  has  dropped  slightly.  This  is  due  to  the  inclusion  of  charged  species.  It  should  be 
pointed  out  that  the  computations  presented  in  these  robustness  studies  were  performed  with 
poor  initial  estimates  in  order  to  test  the  robustness  of  the  Black  Box.  These  extreme  situations 
should  rarely  occur  when  the  Black  Box  is  utilized  in  a  flow  solver,  where  good  estimates  are 
usually  available.  Again,  the  average  number  of  iterations  rarely  exceeds  16. 

6.  EFnCIENCY 

Comparative  timing  runs  were  made  on  a  Cray  XMP  for  the  Mass  Constraint  Technique  and 
Degree  of  Advancement  Technique,  using  the  Ideal  Dissociating  Oxygen  Model,  and  5, 1 1  and 
17-Species  Air  Models.  The  Degree  of  Advancement  Technique  outperformed  the  Mass 
Constraint  Technique:  the  CPU  times  per  iteration  per  computational  cell  were  0.0355, 0.141, 
0.98  and  2.72  msec,  respectively  for  the  Degree  of  Advancement  Technique,  and  were  0.0533, 
0.215, 1.20  and  3.25  msec,  respectively  for  the  Mass  Constraint  Technique.  Thus,  the  increase  in 
the  complexity  of  the  Jacobians  is  less  time  consuming  than  the  computational  savings  afforded 
by  the  reduced  set  of  equations. 

Numerical  experiments  were  performed  in  order  to  test  the  efficiency  of  the  vector  reduction 
implementation.  The  5-Species  Air  Model  was  employed  for  the  computation  of  the  equilibrium 
state  corresponding  to  complete  dissociation  of  oxygen,  with  density  q  =  1.293  10“^  kg/m^ 
and  internal  energy  e  =  6.614  MJ/kg.  The  equilibrium  temperature  is  3996.0  K.  The  initial 
problem  was  composed  of  1000  cells  in  which  the  composition  is  given  by  the  equilibrium  solu¬ 
tion,  and  only  the  initial  estimates  for  temperature  are  perturbed.  Two  cases  were  studied.  The 
first  case  corresponds  to  a  uniform  unconverged  field,  where  all  the  unconverged  cells  start  with 
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the  same  initial  conditions  and  hence  will  have  the  same  convergence  rate.  The  second  case  uti¬ 
lizes  a  variable  unconverged  field,  in  which  the  unconverged  cells  start  at  different  initial  values 
for  temperature,  with  convergence  varying  from  1  to  20  iterations.  For  each  case,  five  computa¬ 
tions  were  performed  for  fields  in  which  the  percentage  of  initially  converged  cells  are  0, 33, 50, 
67, 90  and  99  percent.  All  computations  were  performed  on  a  Cray  XMP,  and  the  cpu  times  were 
recorded. 

The  results  showing  the  vector  reduction  performance  appear  in  Figure  37,  where  the  per¬ 
cent  reduction  in  computational  time  is  given  versus  the  percentage  of  initially  converged  cells. 
The  case  of  the  uniform  unconverged  field  represents  the  baseline  for  all  the  calculations.  The 
initial  overhead  of  the  vector  reduction  implementation,  given  for  the  computation  for  which 
there  are  no  converged  cells,  was  found  to  increase  the  cpu  time  by  about  6  percent.  For  the  vari¬ 
able  unconverged  field,  an  initial  savings  of  30  percent  can  be  realized.  The  break  even  point  for 
the  uniform  unconverged  case  occurs  when  about  9  percent  of  the  field  is  initially  converged. 

For  higher  percentages  of  initially  converged  cells,  vector  reduction  produces  substantial  savings 
in  cpu  time.  This  is  especially  true  for  the  case  of  a  variable  unconverged  field,  which  is  a  more 
realistic  model.  When  used  in  conjunction  with  a  flow  solver,  in  which  a  relatively  high  percent¬ 
age  of  the  computational  cells  are  initially  converged  or  nearly  converged,  a  significant  decrease 
in  the  computational  expense  of  the  Black  Box  should  be  realized.  This  last  statement  will  be 
verified  in  the  results  for  the  flow  solver. 

Numerical  studies  were  also  made  of  the  effects  of  freezing  the  Jacobian  matrix  utilized  for 
the  Newton-Ri^hson  iterations,  with  mixed  results.  Different  strategies  were  employed,  chang¬ 
ing  the  residual  value  at  which  freezing  was  initiated  and  the  number  of  frozen  iterations.  The 
results  showed  that  for  computations  starting  with  a  good  initial  guess  for  the  composition  and 
temperature,  a  fieezing  strategy  could  be  useful.  A  more  detailed  discussion  is  given  in  Cinnella 
and  Cox  [Reference  25]. 


42 


TCfOOOiQ 


Figure  1  Mole  fractions  versus  temperature.  5-Species  Air  Model,  Vibrational  Model, 
1.293  kg/m3. 


T(1000K) 


Figure  2  Mole  fractions  versus  temperature.  S-Species  Air  Model,  Vibrational  Model, 
Q  =  1.293  10-2kg/m3. 
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Figure  3  Mole  factions  versus  temoerature.  5-Species  Air  Model,  Vibrational  Model, 
1.293  10-4  kg/m^. 
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Figure  4  Temperature  versus  internal  energy.  S-Species  Air  Model,  Vibrational  Modri. 
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Figure  5  Isentrqiic  index  versus  temperature.  5-Species  Air  Model,  Vibrational  Model 
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Figure  9  Mole  fractions  versus  temperature.  11-Species  Air  Model,  Vibrational  Model 
Q  =  1.293  10-2  kgAn^. 


Figure  10  Mole  finactions  versus  temperature.  13-Species  Air  Model,  Vibrational  Mod¬ 
el,  p  «  1.293  10-2  kg/m^. 
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TCtOOOK) 


Figure  13  Mole  fractions  versus  tsnperatiire.  17-Species  Air  Model,  Vibrational  Mod¬ 
el,  e  =  1.293  lO"^  kg/in^ 
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Figure  15  Isentropic  index  versus  tonpenUure.  17-Species  Air  ModeU  Vibrational 
Model. 


Figure  16  Speed  of  sound  versus  temperature.  17-Species  Air  Model,  VibnUional 
Model 
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Figure  17  Temperature  versus  inteanal  energy.  All  air  models,  ^tautional  Model, 
0  =  1.293  10-2  kgAn3 


Figure  18  bentropic  index  versus  temperature.  All  air  modds.  Vibrational  Model, 
0-1.293  10-2  kg/in3. 
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Figure  21  Speed  of  sound  versus  tnnperature.  Hydrogea-Air  Combustion  Model, 
Vibrational  Model. 
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Figure  22  Mole  fractions  vereus  tenq)erBture.  Hydrogen-Oxygen  Ccxnbustirai  ModeL 
Vibrational  Model,  q  «  2.8306  kgAn^. 


Figure  24  Tempexature  versus  internal  energy.  Hydrogen-Oxygen  Combustion  Model, 
Vibrational  Model. 
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Figure  25  Isentropic  index  versus  temperature.  Hydrogen-Oxygen  Combustion  Model. 


1.00 


0.5 


1.0 


1.5 


2.0 


Figure  26  Speed  of  sound  versus  tmnperature.  Hydrogen-Oxygen  Combustion  Model. 
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Figure  27  Mole  fractions  versus  temperature.  Hydrocaitx>n  Combustion  Model, 
0  =  3.295  10-2  kg/m^ 
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Figure  28  Tmnperature  versus  internal  energy.  Hydrocarbon  Combustion  Model. 
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Figure  29  Isentropic  index  vmus  temperature.  Hydrcxjarbon  Combustion  Model. 
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Figure  32  Temperature  versus  internal  energy.  Argon  Plasma  Model. 


38 


Figure  33  Isentrq}ic  index  versus  temperature.  Argon  Plasma  Model. 
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Figure  35  Comparison  of  mixture  viscosity  evaluation  methods.  11-Species  Air  Model, 
p  s  1  atm. 


Figure  36  Comparison  of  mixture  thermal  conductivity  evaluation  methods.  1 1-Species 
Air  Model,  p  «  1  atm. 
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1  1  Method  1 

Method  2 

Method  3 

Method  4 

Method  6 

Case  1 

%C 

95.98 

100.00 

97.16 

100.00 

100.00 

%B 

0.07 

0.00 

0.00 

0.00 

0.00 

%W 

2.94 

0.00 

0.00 

0.00 

0.00 

%M 

1.01 

0.00 

2.86 

0.00 

0.00 

Avg.  It. 

12.43 

10.91 

16.66 

12.38 

10.91 

Case  2 

%C 

78.29 

100.00 

96.87 

100.00 

100.00 

%B 

0.38 

0.00 

0.00 

0.00 

0.00 

%W 

16.67 

0.00 

0.00 

0.00 

0.00 

%M 

4.76 

0.00 

3.13 

0.00 

0.00 

Avg.  It. 

12.06 

10.34 

15.08 

11.72 

12.66 

Case  3 

%C 

100.00 

100.00 

100.00 

100.00 

100.00 

%B 

0.00 

0.00 

0.00 

0.00 

0.00 

%W 

0.00 

0.00 

0.00 

0.00 

0.00 

%M 

0.00 

0.00 

0.00 

0.00 

0.00 

Avg.  It 

12.88 

10.60 

16.89 

12.88 

10.60 

Table  2.  Robustness  study.  Ideal  Dissociating  Oxygen  Model,  Mass  Constraint  Technique 
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1  II  Method  1 

Method  2 

Method  3 

Method  4 

Method  6  || 

Case  1 

%C 

95.98 

100.00 

98.67 

100.00 

100.00 

%B 

0.07 

0.00 

0.00 

0.00 

0.00 

%W 

2.94 

0.00 

0.00 

0.00 

0.00 

%M 

1.01 

0.00 

1.43 

0.00 

0.00 

Avg.  It. 

12.43 

10.91 

16.06 

12.38 

10.91 

Case  2 

%C 

78.29 

100.00 

98.67 

68.80 

66.28 

%B 

0.38 

0.00 

0.00 

0.00 

0.00 

%  W 

16.67 

0.00 

0.00 

41.20 

44.72 

%M 

4.76 

0.00 

1.43 

0.00 

0.00  1 

Avg.  It 

12.00 

10.30 

16.31 

8.44 

9.48  1 

Case  3 

%C 

100.00 

100.00 

100.00 

99.83 

99.83 

%B 

0.00 

0.00 

0.00 

0.00 

0.00 

%  W 

0.00 

0.00 

0.00 

0.17 

0.17 

%M 

0.00 

0.00 

0.00 

0.00 

0.00 

Avg.  It. 

12.85 

10.47 

13.49 

12.86 

10.47 

Table  3.  Robustness  study.  Ideal  Dissociating  Oxygen  Model,  Degree  of  Advancement 

Technique 
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Method  1 


Method  2 


Method  5 


Setl 

Set  2 

Setl 

Set  2 

Set  1 

Set  2 

Case  1 

%C 

3.33 

6.00 

98.67 

84.00 

100.00 

100.00 

%B 

46.67 

47.60 

0.00 

0.00 

0.00 

0.00 

%W 

61.00 

46.40 

0.00 

0.20 

0.00 

0.00 

%M 

0.00 

0.00 

1.33 

16.80 

0.00 

0.00 

Avg.  It. 

9.80 

11.60 

13.13 

13.84 

16.01 

16.34 

Case  2 

%C 

86.33 

71.20 

98.00 

90.80 

100.00 

100.00 

%B 

0.00 

10.40 

0.00 

0.00 

0.00 

0.00 

%W 

13.67 

18.40 

0.33 

3.40 

0.00 

0.00 

%M 

0.00 

0.00 

1.67 

6.80 

0.00 

0.00 

Avg.  It. 

17.24 

17.28 

16.63 

16.09 

14.24 

14.69 

Case  3 

%C 

73.67 

62.00 

81.00 

66.20 

100.00 

100.00 

%B 

0.33 

1.60 

0.00 

0.00 

0.00 

0.00 

%W 

26.00 

36.40 

1.33 

1.20 

0.00 

0.00 

%M 

0.00 

0.00 

17.67 

33.60 

0.00 

0.00 

Avg.  It.  1  16.93 

16.72 

13.66 

13.42 

16.60 

16.88 

Table  4.  Robustness  study,  5-Species  Air  Model,  Mass  Constraint  Technique 
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Method  1 

Method  2 

Method  6  || 

Setl 

Set  2 

Set  1 

Set  2 

Set  1 

Set  2  1 

Case  1 

%C 

3.33 

6.00 

99.00 

84.20 

28.67 

34.80 

%B 

46.67 

47.60 

0.00 

0.00 

0.00 

0.00 

%W 

61.00 

46.40 

0.00 

10.80 

71.33 

66.20 

%M 

0.00 

0.00 

1.00 

6.00 

0.00 

0.00 

Avg.  It. 

9.60 

11.17 

12.84 

13.64 

13.08 

14.11 

Case  2 

%C 

86.33 

71.20 

100.00 

91.40 

67.67 

66.20 

%B 

0.00 

10.40 

0.00 

0.00 

0.00 

0.00 

%W 

13.67 

18.40 

0.00 

6.80 

42.33 

44.80 

%M 

0.00 

0.00 

0.00 

1.80 

0.00 

0.00 

Avg.  It. 

16.73 

16.78 

16.76 

16.09 

12.01 

12.28 

Case  3 

%C 

73.67 

62.00 

98.33 

90.00 

49.00 

48.40 

%B 

0.33 

1.60 

0.00 

0.00 

0.00 

0.00 

%W 

26.00 

36.40 

0.00 

1.20 

61.00 

61.60 

%M 

0.00 

0.00 

1.67 

8.80 

0.00 

0.00 

Avg.  It. 

16.36 

16.12 

14.72 

14.60 

10.88 

10.96 

Table  5.  Robustness  study,  5-Species  Air  Model,  Degree  of  AdvaiKrement  Technique 
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Method  1 

Method  2 

Method  6  || 

Setl 

Set  2 

Set  1 

Set  2 

Set  1 

Set  2  1 

Case  1 

%C 

3.00 

3.80 

69.67 

64.40 

97.67 

98.40 

%B 

43.67 

47.00 

0.00 

0.00 

0.00 

0.00 

%W 

46.00 

41.00 

0.33 

0.80 

0.00 

0.00 

%M 

7.33 

8.20 

30.00 

44.80 

2.33 

1.60 

Avg.  It. 

9.89 

14.00 

16.31 

17.76 

18.65 

20.01 

Case  2 

%C 

22.00 

23.40 

36.33 

28.20 

99.33 

99.40 

%B 

19.00 

35.60 

0.00 

0.00 

0.00 

0.00 

%W  1  43.67 

31.00 

0.00 

0.00 

0.00 

0.00 

%M 

16.33 

10.00 

63.67 

71.80 

0.67 

0.40 

Avg.  It. 

15.14 

16.14 

12.72 

16.34 

16.44 

17.06 

Case  3 

%C 

42.67 

27.60 

26.67 

16.20 

100.00 

100.00 

%B 

27.67 

41.00 

0.00 

0.00 

0.00 

0.00 

%W 

20.33 

26.20 

0.00 

0.40 

0.00 

0.00 

%M 

9.33 

6.20 

74.33 

83.40 

0.00 

0.00 

Avg.  It. 

14.90 

16.79 

16.79 

16.43 

16.39 

17.96 

Table  6.  Robustness  study,  11-Species  Air  Model,  Mass  Constraint  Technique 
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Method  1 

Method  2 

Method  5  | 

Setl 

Set  2 

Setl 

Set  2 

Setl 

Set  2 

Casel 

%C 

3.00 

3.60 

67.00 

66.40 

26.67 

30.00 

%B 

44.00 

49.40 

0.00 

0.00 

61.00 

46.40 

%W 

61.00 

44.40 

7.00 

6.00 

22.33 

24.60 

%M 

2.00 

2.60 

26.00 

28.60 

0.00 

0.00 

Avg.  It. 

10.00 

13.44 

17.09 

20.10 

14.24 

18.26 

Case  2 

%C 

22.00 

23.40 

53.33 

37.20 

16.33 

10.40 

%B 

22.00 

37.40 

0.00 

0.00 

29.00 

25.60 

%W 

54.00 

38.00 

0.00 

0.60 

64.67 

64.00 

%M 

2.00 

1.20 

46.67 

62.20 

0.00 

0.00 

Avg.  It. 

14.79 

16.85 

14.59 

16.64 

10.10 

10.66 

Case  3 

%C 

43.00 

27.80 

29.00 

20.00 

10.00 

10.20 

%B 

29.33 

41.00 

0.00 

0.00 

26.00 

23.40 

%W 

26.67 

30.40 

0.00 

0.00 

64.00 

66.40 

%M 

1.00 

0.80 

71.00 

80.00 

0.00 

0.00 

Avg.  It. 

13.90 

14.38 

11.09 

13.95 

9.60 

11.88 

Table  7,  Robustness  study,  11 -Species  Air  Model,  Degree  of  Advancement  Technique 
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Figures?  Vector  reduction  peiformance. 
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SECTION  VI 


MATHEMATICAL  FORMULATION:  FLOW  SOLVER 


The  governing  equations  for  three-dimensional,  compressible,  viscous  fluid  flow  (neglect¬ 
ing  body  forces)  for  a  generalized,  time  dependent,  curvilinear  coordinate  system  have  been 
presented  by  several  authors  [References  26, 27, 28].  In  the  following,  the  density  of  the  fluid 
will  be  denoted  by  Q,  pressure  p,  total  specific  energy  eo ,  thermal  conductivity  ^  viscosity  p 
and  the  Cartesian  velocity  components  will  be  u,  v,  w  for  the  x,  y,  z  coordinate  directions,  re¬ 
spectively.  The  governing  equations,  written  in  strong  consowation  law  form,  read 


,  3F  1  3G  .  3H  _  ,  ^Hy 

^  at)  3^  a|  at)  ’ 


(134) 


where  the  dq)endent  variable  vector  is 


QU 

J  Ov 
QW 


(135) 


the  inviscid  flux  vectors  are 


eu 

QUU  +  Ixp 
QVU  +  lyP 

QWU  +  IzP 

CCoU  +  pdxU  +  lyV  +  ^zW) 

qV 

QuV  +  qxp 
evV  +  r\yp 
QwW  + 

QeoV  +  p(qxU  +  qyv  +  q^w) 


(136) 


(137) 
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H  =  J 


the  viscous  flux  vectors  are 


qW 

QuW  +  ^xp 
QVW  +  ^yp 
QwW  + 

OCoW  +  p(txU  +  CyV  +  CzW) 

0 


Fv  =  J 


‘|2 

+  uT^i  +  vT|2  +  WT^ 


Gv  =  J 


0 

Tni 
Tn2 
Tn3 

+  uT  1  +  vT  2  +  wT, 


'n3 


Hv  =  J 


0 


[(^  +  uTj-i  +  vT,^  +  wT^3 

and  the  Jacobian  of  the  coordinate  transfoimation  is 

ja(x,y,z) 


J  =  det 


In  the  above,  the  curvilinear  coordinates  have  been  defined  as 


(138) 


(139) 


(140) 


(141) 


(142) 
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I  =  I  (  x,y,z,t  )  , 

T1  =  T]  (  x,y,z,t  )  , 

^  ^  (  x,y,z,t  )  , 

t  =  t  . 

and  the  contravariant  velocities  nonnal  to  constant  ^ ,  x]  and  ^  surfaces  are  defined  to  be 

U  =  +  U^x  "b  f"  » 

V  =  tit  +  mix  +  VTly  +  WTlz  , 

W  =  Ct  +  uCx  +  V^y  +  W%z  • 


(143) 


(144) 

(145) 

(146) 


The  heat  flux  and  shear  stress  terms  of  the  viscous  flux  vectors,  given  in  Equations  (139)-(141), 


^  ~  <lx^x  +  —  Txj^x  +  “^xy^y  +  ^xi^z  * 

^  “  ^x^lx  +  %^y  +  <lztlz  »  =  Txy^x  +  "^yy^y  +  'tyzSz  . 


=  qx^x  +  qy^y  +  qzSz 


Tfej  —  txj^x  +  “^yzSy  ^zj^z  , 


~  ‘*'xxqx  f"  “^xy^y  f"  ^xz^z  »  “  txxCx  4”  ‘*'xyCy  4*  ‘>>xaSz  » 

'r,i2  “  "^xy^x  4"  “^yy^y  4"  Xyz^z  »  T^  =  txy^x  4"  T’yyCy  4"  Xy^z  > 


T_3  —  Xxz^x  4"  Xyztly  +  Xa^z  »  “  '*x^x  4"  4’  Xzj^z  * 


(147) 


where 


-  _2„rorfc^u._3u,v3uN  /fc5v.„8v,».3vN 
Xxx  -  3 1*  [  2  (  ^x^  +  qx^  4-  ^x^  )  (  +  r\y^  +  ^y^  ) 

■'yy  =  il*  [  2  ( |y||  +  qy|^  +  )  -  ( Ix|^  4-  qx|j  4-  ^x|f  ) 


/  fc  dW  .  „  dw  .  V  dw  X  T 

-  (  I.  3|- +  <1.  )  ]  . 


»  In  I  2  ( l.^+  r,.^  +  )  -  ( l.^  +  1.#+  ) 


a|  "»n 


a?  '''an 
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Xxy  =  M  (l,|  +  +  Cyf  )  +  (  )  1  . 

=  +  +  +  +  ?xf  )1, 


Xyx  =  M  (  lx|+  1x|^  +  ;x|  )  +  (  ?yf  +  1,|^+  ?yf  )  ]  . 

qx  =  *(?xf +%f+txf  ). 

_  _ if  /■  fc  3T  1  „  3T  I  f-  3T  \ 

qy  -  JG  ( ly-^  ■‘‘  ^yaq  ■'■  ^  ’ 

^  _  ov-  /  fc  3T  .  „  .  J.  ^  . 

qz  *^(^231+  112  ■''  ^*a^  ^  ■ 


(148) 


In  the  above,  the  additional  assumption  that  the  flowfield  is  considered  to  be  a  homogeneous 
mixture  has  been  made,  hence  diffusion  effects  are  neglected.  The  modeling  of  diffusion  would 
require  writing  additional  partial  differential  equations  describing  the  conservation  of  mass  for 
each  of  the  elemental  species  comprising  the  mixture  [Reference  12], 

The  continuity  equation  is  represented  by  the  first  row  in  Equations  (134)  -  (141),  whereas 
rows  two  through  four  correspond  to  the  momentum  equations  and  the  last  row  is  the  energy 
equation.  The  entire  set  of  equations  is  commonly  referred  to  as  the  Navier-Stokes  equations, 
although  incorrectly  so,  since  this  name  originally  defined  only  the  momentum  equations.  How¬ 
ever,  the  author  will  adhere  to  the  accepted  practice.  The  system  of  equations  is  closed  by  the 
“Black  Box”,  which  provides  pressure,  temperature,  viscosity  and  thennal  conductivity  as  func¬ 
tions  of  density  and  internal  energy. 


1.  NONDIMENSIONALIZATION 


The  dimensional  quantities  in  the  above  equations  are  scaled  using  the  following  relations 
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where  the  bar  denotes  the  nondimensional  quantities,  and  L,  Qr<  'I'd  and  Of  are  the  reference 
length,  density,  temperature  and  speed  of  sound,  respectively.  The  reference  viscosity  Pr  and 
thermal  conductivity  used  in  the  previous  relations  are 
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Hf  —  , 


%r  =  lirRr  .  (150) 

where  the  reference  gas  constant  has  been  introduced,  and  is  defined  by  die  relation 

a?  =  RrT,  ,  (151) 

Normally,  viscosity  and  thermal  conductivity  are  scaled  by  introducing  the  Reynolds  and  Prandtl 
numbers,  where  perfect  gas  relations  are  used  to  define  the  latter  [References  26, 27, 28]].  This 
is  not  practical  for  real  gas  flows.  After  drt^ping  the  bar  over  the  scaled  variables,  the  resulting 
set  of  nondimensional  Navier-Stokes  equations  will  q}pear  idratical  to  Equations  (134)  -  (141). 

2.  THIN-LAYER  APPROXIMATION 

Solution  of  the  complete  Navier-Stokes  equations  is  both  computationally  expensive  and 
memory  intensive.  This  is  due  to  the  complexity  of  the  viscous  terms  and  the  large  number  of 
computational  volumes  that  would  be  required  to  resolve  the  mixed  secorxl  derivatives  in  those 
terms.  Boundary-layer  simplifications,  while  reducing  the  complexity  of  the  Navier-Stokes 
equations,  are  not  t^licable  to  flows  where  thoe  is  a  strong  interaction  between  the  inviscid 
and  boundary-layer  regions  [Reference  28]. 

The  concqtt  of  the  thin-layer  i^proximation,  where  the  viscous  terms  comprising  deriva¬ 
tives  in  directions  parallel  to  the  body  are  neglected,  not  only  provides  a  reduction  in  the  com¬ 
plexity  of  the  equations  but  allows  for  the  determination  of  mildly  separated  and  reverse  flow 
regions  [Reference  29].  The  neglected  terms  are  small  in  comparison  with  the  other  viscous 
terms  containing  derivatives  normal  to  the  body  surface,  and  this  makes  the  procedure  justifi¬ 
able.  Moreover,  adequate  resolution  of  these  terms  is  difiicult  to  achieve  unless  grids  which  are 
densely  packed  in  a  direction  parallel  to  the  body  surface  are  employed,  and  this  is  beyond  the 
ctnrent  capability  of  practical  viscous  calculations. 

Neglecting  derivatives  in  the  |-direction,  as  an  example,  the  thin-layer  Navier-Stokes 
(TLNS)  equations  will  read 


,  aF  I  3G  I  3H  _  3Gy  3Hv 
dx  dx\^  d%  aq  ’ 


(152) 


where  the  dqtendent  variable  and  inviscid  flux  vectors  ranain  the  same,  as  givoi  by  Equations 
(135)  -  (138).  However,  the  viscous  flux  vectors  Gy  and  Hy  are  simplified.  They  can  be  written 
in  the  general  form 
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Sv  =  J 


^kl 


^k2 


‘k3 


Cl^  +  uT|ji  +  vTj2  +  wTjj3 


(153) 


where  Sv  denotes  the  viscous  flux  vector  in  the  direction  k.  The  viscous  shear  stress  and  heat 
flux  terms  are  now  given  by 


where 


and 


T|£i  —  Xxxkx  "b  txylCy  "b  « 
Tjj2  =  Txykx  +  "^yyl^y  ‘^y^z  » 
Tjj3  =  "txzkx  +  “tyaky  +  Xx^z  » 


(154) 


Cl|,  —  (Jxl^X  "I"  Qyky  +  qzkz  » 
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] , 

txy  =  H  [  ky|^  +  kx|^  ]  , 
txz  =  I*  [  kz|^  +  kx|^  ]  , 
tyz  =  1*  [  kz|^  +  ky|^  )  , 

qx  =  %  kx|^  , 
qy  =  ite  ky|I  . 

<lz  =  %  kz|^  . 


(155) 


In  die  following,  all  the  viscous  solutions  will  be  obtained  by  numerical  solutions  of  die  TLNS 
equations. 
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SECTION  VII 


NUMERICAL  FORMULATION:  FLOW  SOLVER 
1.  FINITE- VOLUME  DISCRETIZATION 


The  governing  fluid  dynamic  equations,  which  were  presented  in  Section  VI,  can  be  reduced 
to  a  set  of  solvable  algebraic  equations  using  the  finite-volume  technique,  where  the  computa¬ 
tional  domain  is  discretized  by  small,  but  finite,  control  volumes  or  computational  cells. 

Integrating  Equation  (152)  over  the  unit  computational  cell  (A^=Aii=A5=l)  and  noting 
that  the  dependent  variables  become  cell-averaged  values,  the  following  discretized  fonn  of  the 
governing  equations  is  obtained 


^  +  6i(F)  -H  6j(G  -  Gy)  +  6^(H  -  Hv)  =  0  ,  (156) 

where  the  central  difference  operator 

^lO  =  (  ),+i  -  (  )i_i .  (157) 

has  been  introduced,  and  indices  i,  j,  k  correspond  to  the  directions  q,  g,  respectively. 


2.  IMPLICIT  ALGORITHM 


Considering  for  a  moment  only  the  inviscid  fluxes  of  Equation  (156),  a  fully  implicit  for¬ 
mulation  can  be  written  as 

AQ"  =  -  Q"  =  -  At  S"+*  ,  (158) 

where  “  represents  the  time  level  and  At  =  x"'*'  ‘  -  x"  is  the  time  step.  In  the  above,  the 
summation  is  over  all  three  flux  vectors,  where  S  corresponds  to  F,  G  and  H,  and  1  corresponds 
to  i,  j  and  k,  respectively. 

The  inviscid  fluxes  are  nonlinear  functions  of  Q  and  can  be  linearized  in  the  mann«'  used  by 
Beam  arxl  Warming  [Reference  30],  and  Briley  and  McDonald  [Reference  31],  as  follows 

S"+»  =  S"  +  5(Q")  AQ"  ,  (159) 

where 

5(Q")  =  .  (160) 

Applying  this  linearization  to  Equation  (158),  the  Euler-Implicit  difference  formulation  is  ob¬ 
tained 
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(  I  +  At  S  OAQ"  =  -  6,  S"  , 


(161) 


which  is  first-order  accurate  in  time.  The  dot  indicates  that  AQ"  is  included  in  the  operation  of 

6i- 

Similar  considerations  can  be  q)plied  to  the  viscous  flux  vectors,  resulting  in  the  following 
fiilly-implicit  algorithm  for  the  TLNS  equations 


{l  +  At  [6iA  •  +  63(6  -  Bv)  •  +  6k(C  -  Cv)  *  ]}aQ" 
=  -  At[6i(F)"  +  6j(G  -  Gv)"  +  &k(H  "  Hy)"]  , 


(162) 


where  the  Jacobians  of  the  inviscid  and  viscous  fluxes  are 


dG 

c  =  ^ 

dQ  * 

dQ  ' 

(163) 

3Gv 
dQ  ’ 

f'  _ 

3.  FLUX-SPLFT  ALGORITHMS 


The  governing  equations  given  in  Equation  (134)  are  hyperbolic  in  nature  and  consequently 
infonnation  is  propagated  only  in  certain  characteristic  directions.  Numerical  techniques,  where¬ 
by  the  fluxes  are  discretized  in  a  manner  to  allow  information  to  pnpagate  in  the  “correct” 
direction,  are  known  as  upwind  methods.  The  two  most  popular  categories  of  upwind  methods 
are  flux-vector  splitting  and  flux-difference  splitting.  Both  of  these  techniques  involve  the 
discretization  of  the  inviscid  flux  vector  in  one  space  dimension,  and  extensions  to  three  dimen¬ 
sions  are  accomplished  by  considering  three  sq}arate  one-dimensional  problems. 

The  flux-vector  split  scheme  that  is  used  in  this  study  is  of  the  Steger- Wanning  type.  Origi¬ 
nally  developed  for  perfect  gases  [Reference  32],  it  has  been  extended  to  flows  in  chemical  equi- 
lilnum  by  Vinokur  and  Montagnd  [Reference  33].  The  basic  premise  behind  the  scheme  is  that 
flux  vectors  are  split  and  discretized  in  directicHis  corresponding  to  the  sign  of  the  propagating 
wave  speeds. 

The  flux-difference  split  scheme  that  is  employed  is  a  newly  derived  ^proximate  Riemarm 
solvo-  fOT  arbitrary  gas  mixtures.  It  is  based  on  the  perfect  gas  version  developed  by  Roe  [Refer¬ 
ence  34].  Essentially,  this  scheme  involves  the  solution  of  local  Riemarm  problems  at  each  cell 
face,  where  a  discontinuity,  created  by  differences  between  the  left  and  tight  states,  is  assumed  to 
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exist.  In  contrast  to  the  flux-vector  split  scheme,  the  inviscid  flux  vector  is  not  spht  but  recon¬ 
structed  from  the  left  and  right  states. 

The  inviscid  flux  Jacobians  are  analytically  derived  from  the  Steger- Wanning  spht  fluxes  in 
a  manner  similar  to  the  work  of  Belk  [Reference  35],  with  the  exception  that  the  present  for¬ 
mulation  is  based  on  flows  in  chemical  equiliteium.  Both  Steger- Warming  and  Roe  fluxes  can 
be  used  to  determine  the  right  hand  side  of  the  discretized  governing  equations,  where  for  hyper¬ 
sonic  flows  the  former  is  more  robust  and  the  latter  is  more  accurate  [Reference  36]. 

a.  EIGENVALUES  AND  EIGENVECTORS 

Using  the  generic  formulation,  where  S  is  F,  G  or  H  when  k  is  q,  5,  the  inviscid  flux  vec¬ 
tors  can  be  written  as 


S  =  JIVkI 


QPk 

euPt  +  pkx 
Qvp^  +  pky  , 
Qwp^  +  pkz 
ehop,t  -  pkt 


(164) 


where  kx,  ky,  kz,  and  kt  are  the  normalized  metrics  and  are  given  by 

it  = 

ivki  =  ykj  +  kf  +  ki . 

In  the  above,  p,^  is  the  relative  contravariant  velocity  with  respect  to  direction  k,  and  is  defrned 
to  be  the  sum  of  the  absolute  contravariant  velocity  and  the  velocity  due  to  the  time  rate  of 
change  of  the  curvilinear  coordinate  1^  as  follows 

Pij  =  0jj  +  ict  ,  0k  =  ukx  +  vicy  +  wicz  .  (165) 

The  Jacobian  matrices  of  the  inviscid  flux  vectors,  as  well  as  the  eigenvalues  aixl  eigenvec¬ 
tors,  have  been  determined  in  the  manner  of  the  perfect  gas  formulaticms  developed  by  Whitfield  . 
and  Janus  [Reference  37],  with  extensions  to  chemical  equilibrium  made  in  this  work.  The  de¬ 
tails  of  this  determination  are  given  in  Apperxlix  B.  In  summary,  the  eigenvalues  are 
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(166) 


XJ  =  Xj  =  Xj  =  JlVklp^  , 

X{  =  JIVkl(  Pk  +  a  )  , 

Xj  =  JlVkl(  Pk  -  a  )  , 

and  the  right  and  left  eigenvectors  are  given  by  Equations  (B.19)  and  (B.20),  respectively.  In  the 
relations  for  the  eigenvectors,  real  gas  effects  are  accounted  for  in  the  term 

tjj  =  a^  -  h(Y  -  1)  ,  (167) 

where  y  was  defined  in  Equation  (29).  It  is  int^esting  to  point  out  that  for  a  pafect  gas,  y  =  y, 
tp  is  identically  zero,  and  the  eigenvectors  will  identically  reduce  to  their  perfect  gas  counterpart. 

b.  STEGER-WARMING  FLUX-VECTOR  SPLIT 

The  Stego*- Warming  flux-vector  split  scheme  for  a  perfect  gas  utilizes  the  homogeneity 
prc9)erty  to  split  the  fluxes  as  follows 

3  3 

i=l  i=l 

where  corresponds  to  the  positive  and  negative  contributions  of  the  three  distinct  eigenva¬ 
lues  to  the  Jacobian  matrix  S.  liou.  Van  Leer  and  Shuen  [Reference  5]  have  shown  that  the  flux 
vector  for  a  real  gas  no  longer  possesses  the  homogeneous  propoty  and  in  fact  is  composed  of 
homogeneous  and  nonhomogeneous  contributions 

S  =  Sh  +  S'  .  (169) 

They  prqx>sed  the  following  pseudo  splitting  for  the  flux  vector 

S*=Sf+is',  (170) 

where  the  homogoieous  contributions  are  split  according  to  the  “standard”  Steger-Waiming  type 
scheme  and  central  differences  ate  used  for  the  nonhomogeneous  ccmtributions.  Alternatively, 
Vinokur  and  Montagn^  [Reference  33]  have  shown  that  the  ^lit  flux  scheme  derived  using  the 
homogeneity  property  is  just  one  solution  of  an  entire  family  of  one-parameter  flux-vector  blit¬ 
tings.  The  scheme  they  proposed  is  based  on  “gamma”  being  defined  as  the  isentropic  index.  The 
final  result  is  a  generalized  fonnulatitm  of  the  Steger- Wanning  flux-vector  splitting  for  an  aibi- 
trary  gas.  The  latter  method  is  utilized  in  this  study  due  to  the  simplicity  of  the  fonnulaticni. 

'The  blit  fluxes,  for  the  generalized  Steger- Wanning  flux-vecUnr  split  fonnulation  devel¬ 
oped  by  Vinokur  and  Montagnd  [Reference  33],  are  given  by 
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■(^  r  -  i) 


Q 

Q(u  +  akx) 
0(v  +  aky) 
q(w  +  akz) 

Q(ho  +  a0jj) 


Q(u  -  akx) 

1  e(v  -  aky)  j^5(±) 
q(w  -  akz) 
e(ho  -  aOfc) 


where 


xj<±>  =  jivkt&s^-Ss! , 

xj...  =  i,vid»t±^|A±J!! . 

X«±)  =  ^  -  l! 


(171) 


(172) 


In  the  above,  F  is  tiie  isentropic  index  defined  by  Equation  (28).  Similaily  to  what  has  been  seen 
for  the  eigenvectors,  the  above  fonnulas  will  reduce  to  the  perfect  gas  fonnulation,  when  F  s  y- 

As  was  discussed  previously,  the  inviscid  flux  Jacobians  are  constructed  from  the  analytical 
differentiatitm  of  die  split  flux  vectors  in  the  manner  used  by  Belk  [Reference  35].  The  details  of 
the  differentiation  as  well  as  the  resulting  componrats  of  the  Jacobian  matrices  are  given  in  Ap- 
poidix  C.  The  components  of  die  Jacobians  ccxitain  derivatives  of  die  isentropic  index,  which  is 
a  nonlinear  function  of  density  and  tenqierature.  In  die  presoit  study,  these  deiivrfives  are  ne¬ 
glected. 

Discretizatimi  of  the  inviscid  flux  Jaccdiians  at  cell  face  1  +  1/2  is  implemotted  as  follows 
[S  AQj  “  P  M+i  -  Si+J  AQ,  +  Si+j  AQ,+  ,  .  (173) 


As  with  the  split  fluxes,  the  Jacobians  of  die  split  fluxes  can  be  shown  to  reduce  to  dieir  pmfect 
gas  formulation. 


78 


c.  APPROXIMATE  RIEMANN  SOLVER 


The  flux-difference-split  algorithm  used  is  based  on  the  scheme  developed  for  perfect 
gases  by  Roe  [Reference  38].  It  involves  the  solution  of  local  Riemann  problems  at  each  cell  in¬ 
terface,  where  a  left  state  ( •  )i  and  a  right  state  ( *  )r  are  defined  by  extrapolation  of  the  cell-cen¬ 
tered  values  to  the  left  and  right,  respectively. 

Development  of  an  approximate  Riemann  solver  hinges  on  the  determination  of  appropriate 

averages  for  the  eigenvalues  X^,  right  eigenvectors  Tj  and  wave  strengths  such  that  the  follow¬ 
ing  jump  conditions  are  satisfied 

5  5 

=  ISl  =  X °i  f i  ■ 

i=l  i=l 

where  the  cell  interface  states  need  not  be  close  and  thus  |[Q]]  is  arbitrary.  In  the  above,  the  jump 
operator  is  defined  as  the  difference  or  jump  in  values  between  die  right  and  left  states, 
tt(  •  )]1  =  ( '  )r  ~  ( '  )i-  The  averaged  eigenvalues  will  read 

Xj  =  X2  =  X3  =  JlVkip^^  , 

X4  =  JIVkI(  Pk  +  a  )  . 

X5  =  JlVkI(  Pk  -  a  )  , 

where  the  directional  subscript  k  has  been  dropped  for  clarity.  The  averaged  eigenvectors  will 
read 


A  Aa 

aukx 

avkx  +  Qkx 
awicx  —  Qky 

okx^ho  -  a^/(Y  -  l)j  +  e(vkx  -  wky) 

oky 

A  A  •  A  « 

auky  -  eki 

A  A« 

avky 

awky  +  Qkx 

aky|ho  -  aV(Y  ”  1))  +  Q(wkx  -  ukj) 
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1  ‘ 

V  ±  akx 
w  ±  aicz 

A  .A 

ho  ±  a0k 


(176) 


In  the  above,  the  hat  ^  denotes  averaged  terms.  Determination  of  diese  algebraic  averages  is  die 
key  step  in  the  development  of  an  approximate  Riemann  solver.  Abgrall  [Refermce  6]  pointed 
out  that  there  are  miilt4>le  solutions  to  the  system  given  in  Equation  (174)  and  hence  the  aver¬ 
ages  are  not  unique,  as  is  tqiparent  by  the  vaiieQr  of  values  published  in  the  literature  [References 
5-9].  The  approach  taken  in  this  work  is  similar  to  that  of  Abgrall,  however  the  aven^  devel- 
qped  here  are  for  an  arbitrary  chemical  composition.  The  primary  difference  between  the  present 
method  and  others  is  in  the  determination  of  the  pressure  derivative  averages,  where  relatively 
complex  formulas  are  used  by  Vinokur  [Refmence  7]  and  Glaister  [Refoence  8].  Grossman  and 
Walters  [Reference  9]  use  an  approximate  algorithm  which  is  valid  for  near  isentropic  flows. 
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The  present  formulation  utilizes  averages  of  the  ratio  of  enthalpy  and  internal  energy  derivatives 
Y  and  isentropic  index  F,  which  can  be  used  to  constmct  averages  for  the  pressure  derivatives. 

In  the  following,  the  standard  Roe  average,  denoted  by  9i),  will  be  given  as  follows 

=  ^ - 7= - 7= - (178) 

jQi  +  VQr 

and  Q  is  given  by  the  standard  geometric  average,  q  =  JqiQt-  The  averages  are  obtained  through 
the  solution  of  Equation  (174).  The  fuU  derivation  of  the  Roe  averaged  variables  for  the  ^proxi¬ 
mate  Riemann  solver  is  given  in  ^pendix  D.  Hie  averages  for  all  Cartesian  and  contravariant 

velocity  vectors  are  found  to  be  the  standard  Roe  averages;  i.e.  u  =  9ij(u)  and  6^  =  3^5(0,^).  To- 

tal  enthalpy  is  given  in  the  same  manner,  ho  =  3^(ho).  The  speed  of  soimd  a  will  read 


A  2  A  A  A  A 

a  =  FRT  +  (y  -  1)1 


a2  N 

A  Q  A  .  A  A 

ho-|--2!Yft-RT 


i=l 


(179) 


where 


(180) 


In  the  preceding  relations,  the  mixture  gas  constant  R,  temperature  T,  species  internal  eneigy  e^, 
and  species  mass  fractions  Y  j,  are  all  defined  using  the  standard  Roe  average  as  follows 


R  =  %(JR)  ,  f  =  gi(T)  ,  Ci  =  «R>(ei)  ,  Yj  =  ,  (181) 

and  the  barred  terms,  specifically  die  species  specific  l^ats  at  constant  volume  and  mass  fraction 
derivatives,  are  taken  as  integral  averages.  The  simplest  method  of  evaluating  these  integral  av¬ 
erages  is  to  use  the  trt^zoidal  rule,  which  results  in  the  following  arithmetic  averages 
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Cvj  ~  ^  Cvj(Tr)|  , 


(182) 


Again,  the  Roe  averages  given  above  can  be  shown  to  reduce  to  their  perfect  gas  formulations. 

Implementation  of  the  approximate  Riemann  solver  follows  the  methodology  of  Whitfield, 
Janus  and  Simpson  [Reference  39],  where  the  first-order  flux  at  ceU  face  1  +  1/2  is  computed 
fiom  the  following 

^  i«I 

where 

Si  =  S(Qi)  ,  Sr  =  S(Q,+  i)  .  (184) 

Higher-oixier  spatial  accuracy  is  achieved  using  the  flux  interpolation  ^proach,  where  higher- 
order  terms,  involving  the  averaged  eigenvalues,  right  eigenvectors  and  wave  strengths,  given  in 
Equations  (175)  -  (177),  are  added  to  the  first-order  flux  given  above  [Reference  39]. 

4.  MODIFIED  TWO-PASS 

The  discretized  equations.  Equation  (156),  are  advanced  in  time  using  the  modified  two- 
pass  factorization  develqied  by  Whitfield  [Reference  43],  which  is  a  basic  modification  of  the 
standard  two-pass  scheme.  Assuming  that  the  viscous  flux  Jacobians  can  be  split  into  positive 
and  negative  contributions.  Equation  (162)  can  be  written  as 

(I  +  At6,M+  •  +  Ax6jM"  OAQ"  =  -  AtR"  (185) 

where 

6iM+  =  6iA+  +  6j(B+  -  B+)  +  6^(0+  -  C+)  .  (186) 

6,M-  =  6iA-  +  6j(B-  -  By-)  +  iVC"  -  Cy")  ,  (187) 

and  AtR“  represents  the  entire  RHS  of  Equation  (162).  Afte-  expanding  and  manipulating  Equa¬ 
tion  (185),  the  result  can  be  written  in  the  following  factored  form 
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(I  -  D“^M+_i)(I  +  Df ‘M,+  i)AQ"  =  -  ,  (188) 

where 

D,  =  ^  +  M+  +  Mf  .  (189) 

The  above  equation  can  be  solved  in  the  following  two  steps 

(D,  -  M+i)AQ*  =  -  R;*  , 

(Di  +  M,+i)AQ"  =  D,  AQ*  .  (190) 

using  an  efficient  LU  decomposition. 


5.  VISCOUS  FLUXES 

Computation  of  the  viscous  fluxes  follows  the  methodology  of  Chen  [Reference  28],  with 
the  Black  Box  providing  the  necessary  transport  and  thermodynamic  prqserties  in  place  of  the 
perfect  gas  relations.  However,  Chen  [Refererwe  28]  and  others  treat  the  viscous  terms  explicitly, 
since  implicit  treatment  using  a  characteristic-based  splitting  is  not  possible.  In  the  following, 
the  method  used  to  split  the  viscous  fluxes,  in  a  manner  that  allows  implicit  treatment  of  these 
terms  and  preserves  the  efficiency  of  the  existing  LU  algorithm,  will  be  discussed.  Also,  a  brief 
discussion  will  be  given  as  to  the  modeling  of  turbulent  flows. 

a.  VISCOUS  FLUX  DISCRETIZATION 


The  generic  viscous  flux  vector  given  in  Equation  (153)  can  be  rewritten  as 

0 

^  +  Sk.v 


where 


Sv  =  JIAk|2 


11.^  4.  iiv  V* 

+  — k  V* 
•*ak  ^  3^"^ 


ak  ^  ak 


aT 


(191) 


q2  =  u^  +  v^  +  w^  , 
V*  =  k 


(192) 
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A  symbolic  representation  of  the  terms  comprising  the  viscous  flux  at  face  1+1/2  can  be  given 
as  follows 


[Sv(Q)],,.  =  (O), 


L  J.H-i  ■ 


(193) 


In  the  above,  the  flux  is  divided  into  derivative  terms  (  *  )  and  all  others  (  O  ),  where  the  latter 
comprises  transport  properties,  p  and  %,  and  metric  terms  kx. .  The  derivatives  are  discretized 
using  central  differences  as  follows 


All  other  terms  are  evaluated  as  averages 

_(0),  +  (0),+  i 
- 2  ’ 

with  the  excqjtion  of  the  metric  terms,  which  axe  discretized  in  the  following  manner 

(  JlAklkx, 

(Akl(  kx,  *  1 - -  . 

where  J  is  the  Jacobian  of  the  coordinate  transfcsmation. 


(194) 


(195) 


(196) 


b.  VISCOUS  JACOBIANS 

As  was  stated  previously,  the  viscous  fluxes  can  be  linearized  in  a  similar  fashion  as  the  in- 
viscid  fluxes 


S?+^  =  S5  +  5v(Q")  AQ"  , 


where  the  viscous  Jacobians  are  given  as 


Sy(Q")  = 


(197) 


(198) 


IDropping  the  temporal  superscript  and  recalling  the  notation  used  in  the  previous  subsection,  the 
viscous  flux  correction  at  cell  face  1  +  1/2  can  be  written  in  the  following  manner 


3Sy((3) 

aq 


Using  the  same  differencing  techniques  given  above  and  recalling  the  upwind  discretization  of 
the  inviscid  fluxes.  Equation  (199)  can  be  discretized  as  follows 
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(200) 


where 


~3Sv(Q) 

dQ 


AQ 


1+5 


1+i 


Sv  AQj  +  Sy  AQj^  j  , 


(201) 


(202) 


This  is  the  same  type  of  discretization  used  by  Cinnella  [Reference  40]  in  the  develq>ment  of  the 
production  code  GASP  (General  Aerodynamic  Simulation  Program).  The  use  of  Equation  (198) 
allows  for  the  implicit  treatment  of  the  viscous  fluxes  while  maintaining  the  efficiency  of  the  LU 
deconq)osition.  The  discretized  components  of  the  viscous  Jacobians  are  given  in  ^pendix  E. 


c.  TURBULENCE  MODELING 

The  TLNS  equations  are  applicable  to  both  laminar  and  turbulent  flows,  provided  dmt  the 
turbulence  model  is  based  on  eddy  viscosity.  The  total  viscosity  and  thermal  conductivity  are 
given  as  the  sum  of  the  laminar  and  turbulent  contributions,  as  follows 

1*  =  1*1  +  |At  .  (203) 

%  =  +  (204) 

where  subscripts  1  and  t  denote  the  laminar  aixi  turbulent  values,  respectively.  The  lamiruu'  vis¬ 
cosity  and  thermal  coixluctivity  are  provided  by  the  Black  Box. 

The  turbulent  viscosity  is  determined  using  the  Baldwin  and  Lomax  model  [Reference  41]. 
A  very  good  description  of  this  model  and  its  application  is  given  by  Chen  [Reference  28]  and 
the  details  are  not  repeated  here.  The  turbulent  thermal  coixluctivity  is  computed  by 

as,  =  ^  .  (205) 

where  Cp  is  the  frozen  specific  heat  at  constant  pressure  defined  in  Equation  (19)  and  Pit  is  the 
turbulent  Praixltl  number.  Since  there  are  few  turbulent  studies  for  high  enthalpy  flows  aixl  no 
data  as  to  what  is  a  practici^l  value  for  Prt,  the  staixlaid  value  of  Prt  =  0.90  will  be  used  in  the 
calculations. 
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6.  BOUNDARY  CONDITIONS 


Explicit  boundary  conditions  are  applied  throughout  using  a  first  order  extrapolation  to  a 
layer  of  phantom  cells  surrounding  the  computational  domain.  All  inflow  and  outflow  bound¬ 
aries  incorporate  the  characteristic  variable  boundary  conditions,  CVBC’s,  as  developed  by  Janus 
[Reference  42],  which  allow  flow  information  to  be  prc^agated  in  the  direction  of  the  character¬ 
istic  wave  speeds.  CVBC’s  are  also  used  at  impermeable  surfaces  for  inviscid  computations.  No¬ 
slip  boundary  conditions  are  used  at  impermeable  surfaces  for  viscous  computations  along  with  a 
specified  pressure  gradient  and  either  adiabatic  wall  or  specitied  wall  temperature. 

Application  of  the  boundary  conditions  for  a  real  gas  are  the  same  as  for  a  perfect  gas,  with 
the  exception  that  the  thennodynamic  properties  are  obtained  by  means  of  the  Black  Box,  as  was 
stated  previously  in  Section  n.  The  known  thermodynamic  properties  resulting  hrom  the  enforce¬ 
ment  of  the  CVBC  subsonic  inflow,  outflow  or  impermeable  surface  boundary  condition  are 
pressure  and  density  and  therefore  the  complete  thermodynamic  system,  CTS,  is  detined  by  us¬ 
ing  the  “Black  Box”  as  a  function  of  pressure  and  density,  CTS  =  f(p,  Q).  For  the  no  slip  bound¬ 
ary  condition  with  specified  temperature,  pressure  and  temperature  are  the  known  properties  and 
hence,  CTS  =  f(p,T). 
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SECTION  vin 


FLOW  SOLVER  RESULTS 

Numerical  results  have  been  obtained  using  the  methods  proposed  in  the  previous  sections. 
As  a  means  of  validation,  four  test  cases,  two  inviscid  and  two  viscous,  were  examined.  The  in- 
viscid  computation  involved  hypersonic  flow  about  a  blunt  9®  half-angle  cone  with  a  nose  radius 
of  2.5  in.  and  the  numerical  simulation  of  the  Space  Shuttle  Main  Engine  (SSME).  The  two  vis¬ 
cous  test  cases  include  a  standard  flat  plate  computation,  for  code  validation,  and  a  hypersonic 
inlet.  Also,  a  numerical  study  is  made  of  the  efficiency  of  the  chemical  equilibrium  solver. 

All  of  the  solutions  have  been  obtained  using  second  order  spatial  accuracy,  advancing  die 
solution  to  its  steady  state  by  means  of  the  Euler  implicit  scheme  described  in  Section  VII.  The 
approximate  Riemann  solver  of  the  Roe  type,  developed  in  Section  Vn,  is  used  for  all  calcula¬ 
tions  and  higher  order  accuracy  is  attained  using  the  Van  Leer  limiter. 

Characteristic  variable  boundary  conditions  are  used  wherever  i^iphcable,  including  die  im¬ 
permeable  surface  for  inviscid  computations.  No-slip,  adiabatic-wall,  with  zero  pressure  gradi¬ 
ent  boundary  conditions  are  used  for  the  viscous  calculations. 

The  simple  S-Species  Air  Model  is  employed  as  the  air  chemistry  model,  along  with  the 
Consistent  Kg.  The  thennodynamic  model  used  is  the  Vibrational  Model.  Newton  Limiter  and 
Relative  Limiter  are  used  for  corrections  of  the  species  densities  and  temp^ature  updates,  re¬ 
spectively.  Transport  properties  are  provided  by  the  combination  of  Sudieiland  Model  for  the 
species  values  and  Wilke’s  Rule  for  the  mixture  values. 

1.  INVISCID  RESULTS 

a.  BLUNT  CONE 

The  geometry  of  the  blimt  cone  and  the  71  x  26  grid  used  for  the  inviscid  computations  can 
be  seen  in  Figure  38.  The  fieestream  conditions,  corresponding  to  an  altitude  of  10  km,  were; 
pressure  poB  =  26.5  kPa,  density  =  0.414  kg/m^  and  temperature  T.  =  223  K.  Com¬ 
putations  were  made  at  Mach  number  M«  =  10  and  zero  angle  of  attack  using  both  the  perfect 
gas  model  and  the  5-Species  Air  Model. 

The  effects  of  a  “real  gas’’  are  immediately  apparent  in  Figure  39,  where  the  temperature 
profiles  along  the  stagnation  streamline  are  plotted  for  both  reactive  and  perfect  gas  computa¬ 
tions.  Shock  standoff  distance  and  temperature  levels  are  strongly  affected  by  the  thermochemi¬ 
cal  model  chosen.  The  reduced  temperatures  for  the  chemical  equihbiium  solution  are  due  to  the 


87 


conversion  of  some  of  the  kinetic  energy  to  dissociation  energy,  where  for  a  perfect  gas  all  of  the 
kinetic  energy  would  go  to  thermal  energy.  Higher  densities  in  the  shock  layer  for  the  real  gas 
solution,  characteristic  of  reactive  flows  [Reference  44],  are  seen  in  Figure  40  where  the  density 
profiles  along  the  stagnation  streamline  are  plotted.  Figure  41  depicts  the  temperature  profiles 
along  the  body  surface,  where  large  temperature  differences  between  the  perfect  gas  and  real  gas 
solutions  exist  at  the  nose,  but  the  solutions  are  in  good  agreement  in  the  cooler  regions,  starting 
at  the  shouldo*.  The  previous  results  conelate  nicely  with  the  plot  of  the  mixture  composition 
along  the  body  surface  seen  in  Figure  42.  The  plot  indicates  strong  oxygen  dissociation  at  the 
nose,  but  as  the  flow  moves  along  the  body  and  cools,  recombination  occurs.  Figure  43  depicts 
the  pressure  coefficient  along  the  body  surface.  It  can  be  seen  that  pressure  is  not  strongly  af¬ 
fected  by  real  gas  effects,  due  to  its  “mechanical”  nature  [Reference  45]. 

b.  SSME  NOZZLE 

The  geometry  of  the  SSME  nozzle  and  the  88  x  31  grid  used  for  the  inviscid  computations 
are  shown  in  Figure  44,  Combustion  chamber  conditions  correspond  to  100%  power  at  sea  level 
and  were:  mixture  ratio  of  6.0,  Mach  number  Me  =  0.2,  temperature  Tc  =  3639.0  K,  and  pres¬ 
sure  Pc  =  20.24  MPa.  Computations  were  made  using  both  the  perfect  gas  model,  with  y  =  118, 
and  Hydrogen-Oxygen  Combustion  Model. 

Density  profiles  along  the  centerline  and  wall  of  the  nozzle  are  given  in  Figure  45.  In  tl^ 
following,  the  axial  distances  shown  in  the  figures  will  be  referenced  to  the  throat  position.  The 
results  show  that  there  is  little  difference  between  the  perfect  gas  and  real  gas  solutions.  A  great¬ 
er  disparity  is  seen  for  the  temperature  profiles  in  Figure  46.  Higher  temperatures  for  the  real 
gas  solutions  occur  at  a  point  just  past  the  throat.  The  reason  for  this  becomes  apparent  when 
looking  at  Figure  47,  where  the  composition  along  the  wall  is  given:  as  the  flow  expands  on  the 
downstream  side  of  the  throat,  recombination  reactions  occur,  which  are  exothermic  in  nature. 

Similar  computations  were  made  by  Wang  and  Chen  [Reference  13],  using  a  pressure-based 
flow  solver  and  equilibrium  chemistry.  Comparisons  of  the  real  gas  solutions  obtained  in  the 
present  study  with  the  results  obtained  by  Wang  and  Chen  [Refereix»  13]  are  made  in  Figure  48 
and  Figure  49,  where  plots  of  Mach  number  profiles  and  pressure  decays  are  given,  respectively. 
The  figures  show  excellent  agreement  between  the  two  sets  of  results,  with  the  exception  being 
the  throat  region.  The  present  study  appears  to  do  a  better  job  of  ct^turing  the  nozzle  shock  at 
the  wall.  The  specific  impulse  of  the  SSME  nozzle  was  computed  to  be  460.0  sec,  which  is  also 
in  excellent  agreement  with  the  value  of  460.4  sec  predicted  by  Wang  and  Chen  [Reference  13]. 
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2.  VISCOUS  RESULTS 


a.  FLAT  PLATE 

In  order  to  validate  the  viscous  c^ability  of  the  code,  solutions  for  the  laminar  and  turbu¬ 
lent  flow  over  a  flat  plate  were  obtained.  The  freestream  conditions  were:  Mach  number 
M  ao  =  0.5,  temperature  T«  =  300  K,  and  Reynolds  numbers  per  unit  length.  Re  /L,  of  1.972 
10^  and  1.972  10®,  for  the  laminar  and  turbulent  computations,  respectively.  The  laminar  results 
for  both  the  perfect  gas  and  5-Species  Air  Model,  presented  in  Figtire  50,  nicely  reproduce  the 
Blasius  solution.  Similarly,  the  turbulent  results  using  the  5-Species  Air  Model,  as  given  in  Fig¬ 
ure  51,  are  seen  to  be  in  excellent  agreement  with  the  law  of  the  Wall  results. 

b.  INLET 

The  hypersonic  inlet  is  composed  of  a  flat-surfaced  cowl  and  a  10°  ramp,  with  an  inlet 
height  D  =  2  m  and  exit  height  of  0.5  D.  The  overall  length  is  4.25  D  and  tte  cowl  length  is  0.5 
D.  Conditions  at  the  inlet  are:  Mach  number  Moo  =  10,  pressure  poo  =*  1.015  kPa,  density 
e  00  =  1.172  10“^  kg/m^,  temperature  T*  =*  300  K,  and  Reynolds  number  based  on  the  inlet 
height  Re^  =  4.4  10®.  Figure  52  depicts  the  inlet  geometry  and  111  x  61  grid  used  for  the  vis¬ 
cous  calculations.  Computations  were  made  using  the  perfect  gas  model  and  die  5-Species  Air 
Model.  Grid  packing  at  the  wall  resulted  in  y  values  less  than  1,  which  is  considered  more  than 
adequate  to  define  the  viscous  sub-layer. 

Figure  53  depicts  the  temperature  profiles  along  the  tq)  and  bottom  walls  of  the  inlet.  Two 
important  points  need  to  be  made  about  this  plot  The  first  is  the  gross  overproduction  of  the 
temperature  values  by  the  perfect  gas  model.  The  second  is  the  displaced  location  of  the  incident 
shock  on  the  top  wall;  shock  locaticxi  for  the  real  gas  solution  is  slightly  downstream  of  that  for 
the  perfect  gas  solution.  The  surface  pressure  distributions  for  die  top  and  bottom  walls  are  given 
in  Figure  54.  As  was  seen  previously,  pressure  is  not  strongly  affected  by  the  inclusion  of  real 
gas  effects.  Similar  results  are  seen  in  Figure  55,  where  the  coefficient  of  skin  friction  is  plotted. 
The  primary  difference  in  both  plots  is  die  displacement  of  the  incident  shock.  The  mixture  com¬ 
position  along  the  bottom  wall  is  given  in  Figure  56,  where  it  can  be  seen  that  oxygen  is  almost 
entirely  dissociated.  Of  particular  interest  is  die  slight  decrease  in  the  mole  percentage  of  the  mo¬ 
natomic  species,  which  occurs  before  the  ramp  irKluced  oblique  shock.  This  correlates  nicely 
with  the  slight  drop  in  temperature  seen  in  Figure  53.  Opposite  results  occur  for  the  expansion. 
Figure  57  depicts  the  Mach  number,  temperature  and  pressure  profiles  along  the  centerline.  As 
can  be  seen,  there  is  very  htde  difference  between  the  solutions  for  the  two  thermochemical 
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models,  other  than  the  reduced  shock  angle.  The  real  gas  effects  are  local  to  the  boundary  layer 
and  have  the  greatest  impact  on  the  surface  temperatures. 

3.  PERFORMANCE 

In  order  to  deteimine  the  efficiency  of  the  chemical  equilibrium  flow  soIvct,  a  numerical 
study  was  made  involAong  inviscid  computations  for  the  blunt  cone  geometry.  CPU  times  and 
iteration  counts  obtained  using  the  5-Species  Air  Model  were  compared  to  those  obtaiiKd  using 
the  perfect  gas  model.  Several  techniques  for  enhancing  the  efficiency  of  the  solver  were  imple¬ 
mented  and  tested.  All  computations  were  made  on  a  Silicon  Gn^hics  IRIS  4D-340VGX  ma¬ 
chine. 

The  techniques  used  to  enhance  the  efficieiKy  include  the  use  of  vector  reduction  and  chem¬ 
istry  freezing,  both  of  which  have  been  described  previously  in  Section  IV.  Five  test  cases  were 
tried  using  the  5-Species  Air  Model  and  compared  to  the  perfect  gas  solution.  Case  1  represents 
the  baseline  computation,  where  no  efficiency  techniques  are  implemented.  Vector  reduction  was 
implonented  in  Cases  2-5,  and  various  freezing  strategies  were  employed  for  Cases  3-5.  The 
convergence  strategy  used  for  both  perfect  gas  and  real  gas  runs,  as  well  as  the  efficiency  strate¬ 
gies  employed  are  summarized  in  Table  8. 

The  results  for  Cases  1-3  and  the  perfect  gas  run  are  given  in  Figure  58,  where  the  normal¬ 
ized  residuals  versus  CPU  time  are  plotted.  Each  symbol  on  the  plots  corresponds  to  100  cycles. 
The  perfect  gas  computation  achieves  a  residual  reductitm  of  almost  five  orders  of  magnitude  in 
a  little  less  than  one  CPU  hour.  The  baseline  real  gas  computation.  Case  1,  requires  20  CPU 
hours  to  attain  a  comparable  residual  reduction.  It  should  be  pointed  out  that  the  first  200  cycles 
accounts  for  nearly  30  percent  of  the  total  CPU  time.  Using  vector  reduction,  a  reduction  in  CPU 
time  of  almost  one  half  is  realized,  as  can  be  seen  for  Case  2.  A  more  significant  reduction  is  at¬ 
tained  when  the  first  freezing  strategy.  Case  3,  is  employed,  resulting  in  a  CPU  time  savings  of 
about  85  percent  over  the  baseline  computation. 

Slight  improvements  are  still  obtainable,  as  is  seen  in  Figure  59,  in  which  the  results  for 
Cases  3-5  are  ctxnpared  to  the  perfect  gas  result  The  total  CPU  time  required  for  Case  4  is  just 
under  three  hours.  However,  CPU  requirements  appear  to  be  little  affected  by  the  particular 
freezing  strategy  employed,  provided  that  a  freezing  strategy  is  employed.  The  results  for  Case 
5,  where  a  marginal  improvement  over  Case  4  can  be  seen  for  the  first  200  cycles  but  no  overall 
gain  is  achieved,  ctmfiim  this  conclusion.  Again  it  should  be  pointed  out  that  the  begiiming  of 
the  computation,  when  variations  in  flow  conditions  are  significant,  requires  a  large  pocentage 
of  die  overall  CPU  time. 
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Figure  38  Blunt  cone  geometry  and  inviscid  grid  (71x26) 


Figure  39  Temperature  profiles  along  the  stagnation  streamline.  Blunt  cone,  invisicid. 


Figure  42  Mole  fractions  versus  distance  along  the  body  surface.  Blunt  cone,  inviscid 
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Figure  44  SSME  nozzle  geometry  and  inviscid  grid  (88x31) 
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Figure  48  Comparison  of  Mach  number  distributions.  SSME  nozzle. 
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Figure  50  Laminar  velocity  profiles  for  the  flat  plate. 
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Figure  56  Composition  along  the  bottom  wall.  Inlet,  viscous. 


Figure  57  Mach  number,  pressure  and  temperature  distributions  alcmg  the  coiteriine.  Inlet,  viscous. 
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Cycle 

CFL# 

Order 

Vector  Reduction  /  Niunber  of  Cycles  Frozen 

Case  #1 

Case  #2 

Case#3 

Case  #4 

Case  #6 

0-100 

0.6 

1st 

N/0 

Y/0 

Y/6 

Y/10 

Y/20 

100-200 

1.0 

2nd 

N/0 

Y/0 

Y/10 

Y/20 

Y/20 

200-300 

2.0 

2nd 

N/0 

Y/0 

Y/10 

Y/20 

Y/60 

300^00 

4.0 

2nd 

N/0 

Y/0 

Y/26 

Y/60 

Y/60 

400-1000 

4.0 

2nd 

N/0 

Y/0 

Y/60 

Y/60 

Y/60 

Table  8  Convergence  and  efficiency  strategies. 


Figure  58  Comparison  of  residuals  versus  CPU  time.  Cases  1-3. 
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Figure  59  Comparison  of  residuals  versus  CPU  time.  Cases  3-5. 
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SECTION  IX 


SUMMARY  AND  CONCLUSIONS 

The  goal  of  the  present  study  has  been  the  development  of  an  efficient  solver  for  flows  in 
local  chemical  equilibrium.  The  work  has  been  carried  out  in  two  phases:  the  development  of  a 
“Black  Box”  solver  for  the  determination  of  the  equilibrium  composition  and  thermodynamic 
properties;  and  the  incorporation  of  this  Black  Box  into  a  perfect  gas  flow  solver  that  has  been 
modified  to  include  real  gas  effects.  The  end  result  is  a  numerical  tool,  whereby  the  solutions  for 
high-temperature,  high-speed  flow  fields  can  be  efficiently  obtained  without  the  enormous  com¬ 
putational  cost  associated  with  finite-rate  chemistry. 

A  Black  Box  chemical  equilibrium  solver  was  developed  for  homogeneous  mixtures  of  ther¬ 
mally  perfect  gases.  Several  practical  chemistry  models,  describing  air  mixtures,  combustion 
mixtures  and  a  plasma,  were  implemented.  The  solver  was  written  in  a  generic  fashion,  whereby 
the  inclusion  of  additional  chemistry  models  of  interest  is  a  relatively  simple  process.  Two  solu¬ 
tion  methods  for  the  determination  of  the  equilibrium  composition  were  derived  and  implement¬ 
ed,  the  Mass  Constraint  Technique  and  the  Degree  of  Advancement  Technique.  A  third  algo¬ 
rithm  was  explored  briefly,  but  was  found  to  be  inadequate  for  the  problem  at  harxi,  although 
further  testing  may  be  indicated.  The  Black  Box  provides  the  choice  between  two  thermodynam¬ 
ic  models.  Vibrational  Model  and  Curveflt  Model,  and  two  methods  for  the  evaluation  of  the 
equilibrium  constant,  Curvefit  Kc  and  Consistent  IQ.  For  viscous  computations,  two  procedures 
for  the  evaluation  of  the  species  transport  properties  were  implemented,  Sutheiiand  Model  and 
Gupta  Model,  aixl  mixture  values  can  be  obtained  with  either  Wilke  Rule  or  Gupta  Rule. 

The  Black  Box  was  thoroughly  tested  fen*  both  efficiency  and  robustness.  Although  short¬ 
comings  were  found  in  the  basic  algorithms,  several  efficieiKy  saving  and  robustness  enhancing 
techniques  were  implemented.  The  end  result  is  a  reliable  solver  for  the  determination  of  the 
composition  and  thermodynamic  properties  for  an  equilibrium  gas  mixture.  The  chemistry  mod¬ 
els  that  have  been  implemented  were  also  tested  on  various  conditions.  The  results  show  that  the 
Black  Box  accurately  predicts  both  composition  and  properties  for  a  wide  range  of  chemistry 
models,  arrd  illustrate  the  limitations  of  the  perfect  gas  assumptioa 

The  inemporation  of  “real  gas”  effects  into  an  existing  perfect  gas  flow  solver  posed  an  ex¬ 
citing  challenge.  Several  modifications  to  the  existing  code  were  made,  including  the  develop¬ 
ment  of  a  new  approximate  Riemann  solver  for  an  arbitrary  real  gas.  Real  gas  versions  for  the 
Steger- Warming  flux-vector  split  scheme,  split-flux  Jacobians,  and  Roe-type  flux-difference 
split  scheme  were  derived  and  implemented.  The  real  gas  formulations  were  develqred  in  a 
manner  whereby  they  reduce  to  their  perfect  gas  versions  in  a  straightforward  fashion. 
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Hypersonic  flows,  both  inviscid  and  viscous,  for  a  few  test  cases  were  investigated,  as  well 
as  a  validation  computation  using  the  ubiquitous  flat  plate.  The  results  were  encouraging,  and 
demonstrated  the  fallacy  of  relying  on  the  perfect  gas  assumption  for  high  Mach  number  flows. 
The  results  for  a  hypersonic  inlet  also  showed  the  importance  of  viscous  effects  in  the  study  of 
real  gas  flows.  Finally,  an  efficiency  study  was  made  of  the  chemical  equilibrium  flow  solver. 
Techniques  implemented  to  enhance  efficiency  worked  very  well.  The  results  indicated  that 
while  the  real  gas  computations  are  more  expensive  than  comparable  perfect  gas  computations, 
they  are  within  reason  and  much  more  so  than  finite-rate  computations. 

Further  work  is  necessary  in  several  areas.  Increasing  the  robustness  of  the  Degree  of  Ad¬ 
vancement  Technique  could  also  increase  the  efficiency  of  the  flow  solver,  since  this  technique 
was  found  to  be  the  more  efficient  of  the  “Black  Box”  algorithms.  Further  investigation  of  the 
chemical  equilibrium  solver  developed  by  Mein^es  and  Morgan  [Reference  4],  might  also  prove 
profitable.  An  extended  study  of  the  effects  of  the  number  of  species  on  flow  solver  efficiency  is 
also  warranted.  In  the  present  study,  the  “gamma”  derivatives  appearing  in  the  split  flux  Jaco- 
bians  were  assumed  to  be  zero,  and  the  effect  of  this  need  to  be  ascertained.  Numerical  studies 
also  need  to  be  made  of  the  transport  property  calculation  procedures.  The  effects  of  using  an 
improved  initial  guess  for  the  Black  Box  computations,  discussed  in  Section  IV,  also  need  to  be 
investigated.  Finally,  multiblock  capability  needs  to  be  added  to  the  code,  providing  the  ability  to 
obtain  solutions  for  more  challenging  geometries,  such  as  the  national  aerospace  plane  (NASP), 
space  shuttle  and  hypersonic  vehicles  with  stores. 
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APPENDIX  A 

SPEED  OF  SOUND  AND  PRESSURE  DERIVATIVES 
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1.  SPEED  OF  SOUND 


The  speed  of  sound  is  defined  by  the  relation 

=  (I).  •  (A.1) 

Writing  the  combined  1st  and  2nd  Laws  of  Thmmodynamics,  for  a  gas  in  chemical  equihbrium 

Tds  =  de  -  ^dQ  ,  (A.2) 

the  partial  derivative  of  internal  energy  with  respect  to  density  at  constant  entropy  can  be  ob¬ 
tained  as 


(A.3) 


Expanding  the  partial  derivative  in  Equation  (A.1)  and  using  the  result  firom  Equation  (A.3),  the 
equation  for  the  speed  of  sound  can  be  written  in  the  following  form 


[dQf^  Q^\del^  • 


(A.4) 


The  thermal  and  caloric  equations  of  state  for  an  equilibrium  mixture  of  thermally  perfect  gases 
can  be  written  as 


N 

P  =  qRT=  ^QiRjT, 
i=l 

N  ^  r 

i=i  i-1  [ 

The  above  equations  can  be  cast  in  differential  form,  as  follows 

N 

dp  =  QRdT  +  qT  ^  RjdT  -1-  RTdg  , 

N 

de  =  CydT  +  ^ejdYj  , 

i=l 

where  the  differential  of  the  mass  fractions  is  given  by 


(A.5) 

(A.6) 

(A.7) 

(A-8) 

(A.9) 
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Substituting  Equations  (A.8)  and  (A.9)  into  Equation  (A.7),  the  pressure  differential  can  be  writ¬ 


ten  as 


N 

qR  +  qT  ^ 


dp  = 


Zj^i\erT) 


+  RT  + 


-Hm 

^  /avA] 


(A.  10) 


Now  the  partial  derivatives  of  pressure  with  respect  to  internal  energy  at  constant  density  and 
with  respect  to  density  at  constant  internal  energy  can  be  easily  obtained  as 


[del 

\  /  O 


R  +  T 


N  ,  , 

i=l _ 

N  ,  , 

Z-  ®>v  ^  / 


(A.11) 


/ao\  «  /aY\  " 

ft .  •  «.  S  -  t  r  ■■ 


(A.  12) 


Substituting  these  results  into  Equation  (A.4)  and  after  some  algebra,  the  equation  for  the 
speed  of  sound  can  be  written  in  the  more  familiar  foim 


-  r| . 


(A.13) 


where  the  isentropic  index  is  defined  as 


107 


R  +  T  X  Ri(^) 

_ i=J^ _ 

N  ,  . 

Cv  +  ^  ar  j 

i=l 

R  +  T  y  Rj® 


(A.  14) 


Defining  y  to  be  the  ratio  of  the  enthalpy  and  internal  energy  derivatives  at  constant  density,  as 
follows 


^=1  + 


(^)  r.t|r.(^) 

-r-T^  =  1  + - — - 

\?FL  ;;  . 

V  Cv  -I- 


(A.15) 


and  substituting  this  definition  into  Equation  (A.  14),  the  isentropic  index  can  be  rewritten  as 


r  =  Y +  ^ 

RT 


i(^)  [V  -  (7  -  l)ej  ■ 


(A.  16) 


Similaily,  the  pressure  derivatives  defined  in  Equations  (A.11)  and  (A.12),  can  be  rewritten  as 


\dej 

\  /q 

/'ip\  ^ 

W). 


e(Y  - 1) , 


RT(1  +  r  -  Y) 


(A.17) 


(A.18) 


2.  PRESSURE  DERIVATIVES 

A  key  part  in  Uk  develqrment  of  a  flow  solver  for  real  gases  is  the  evaluation  of  the  pres¬ 
sure  derivatives,  specifically  the  partial  derivatives  of  pressure  with  respect  to  the  dqjendent 
variable  vector,  Q,  as  follows 


/iP\  =  r lJL.\  / 

\dQ)  [\ae;’  \aQu;’  \aev;’  \aQw/’  \aeeo/ J  • 


(A.19) 


The  pressure  differential  can  be  written  as 


108 


(A.20) 


dp 


de 

Q 


=  +  (^)deu  +  (^)dev  +  (^)dew  +  (5^)00*^ , 

and  internal  energy  is  given  by 


efo  _  (Qu)^  +  (Qv)^  +  (QW)^ 

e  20^ 


(A.21) 


Differentiating  Equation  (A.21),  and  substituting  the  result  into  Equation  (A.20),  the  following 
relations  for  the  pressure  derivatives  are  found 


/  3p\  _  /  3p\ 

\aej  “  \aQj^ 

(ae)^  »(2  '')  +  “*• 

(A.22) 

/  \  ^  _  u/^ 

\dQ\ij  e\ae 

]  =  -  u(Y  -  1)  , 

'  Q 

(A.23) 

(Jr) 

[dQVj 

\  ^  v/5p 
r  o\ae 

1 

> 

1 

II 

(A.24) 

(Jr) 

1  _  W  /5p 

f  e\ae 

II 

1 

1 

(A.25) 

fJE_) 

\dQeo) 

-  1 
eUe/ 

\  /Q 

=  Y-  1  , 

(A.26) 

where  q^  =  u^  +  v^  -t-  w^. 


APPENDIX  B 

EIGENVALUES  AND  EIGENVECTORS  OF  THE  EMVISCID  FLUX  JACOBIANS 
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The  following  development  parallels  the  work  done  by  Janus  [Reference  42]  and  therefore 
many  steps  will  be  intentionally  left  out.  The  primary  focus  is  the  extension  of  the  previous  work 
to  real  gases. 


1 .  CONSERVATIVE  VARIABLE  JACOBIAN  MATRIX 

The  inviscid  flux  Jacobians  are  given  by 

c  -  iS 
^  dQ  ’ 

where  the  dependent  variable  vector  and  generic  flux  vector  are 


Q=  J 


e  ■ 

Ou 

QV 

,  S  =  J 

0W 

QCo 

ePk 

QUpk  +  pkx 

QvPk  +  Pky 
Qwp^  +  pkz 
PkCeCo  +  p)  -  pkt 


(B.l) 


(B.2) 


After  differentiating,  the  components  of  the  inviscid  Jacobian  matrix  are  given  by  rows  as 

“  kx  ,  S^j2  “  ky  ,  S|3  —  kz  ,  =  kj  ,  Sj5  =  0  ,  (B-3) 


S21  - 


S22  “ 


S24  - 


^31  ” 


/ 3p\  _  / 3p\ 

\ /  \ 3e I 

-S(S)J 

/£p\  k 
[del 

( 3p\  _  /3p\ 

1 3o /  1 3e I 

/£p\ 

\dej 


q*' 

" 

\ 

®  2 

I 

C 

'Q 

kx  -  uOu  , 


ukx  +  Pk  , 


=  - 


[dej 

\  /  O 


+  ukv 


_  w/£P 
Q 


+  uk  S  =  —  —( — ^  k 
T  UKz  ,  O25  -  Q  I  ae  /  ’ 

\  /ft 


e  -ST 

® _ 2 


ky  -  V0,,  , 


_  il/£E 
0 


^32  “  O  I  ’  ^33 


^34  = 


_w/£pW  .  5  =_i/£E\k 

eUei ^35  Qlael'^y’ 

\  /  o  '  '  D 


Vky  +  Pjj  , 


(BJ) 
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The  main  difference  between  the  present  fonnulation  and  the  perfect  gas  fonnulation  is  die  eval- 
uaticm  of  the  pressure  derivatives.  Substituting  die  relations  for  the  pressure  derivatives  devel¬ 
oped  in  ^ipendix  A,  the  second  through  fifth  row  (Equations  (B.4)-(B.7))  of  the  Jacobian  ma¬ 
trix  can  now  be  rewritten  as 

S21  “  (^  ''p)kx  “  , 

S^22  =  (2  -  Y)ukx  +  Pk  .  S23  =  -  v{7  -  l)kx  +  uky  ,  (B.8) 

^24  =  -  w(y  -  l)kx  +  uk*  ,  S23  =  (y  -  l)kx  . 

^31  “  (4^  “  ^k  » 

^32  =  “  u(y  -  l)ky  +  vkx  ,  ^33  =  (2  -  Y)vky  +  Pk  ,  (B.9) 

^34=  -  w(y  -  l)ky  +  vki  ,  ^35  =  (y  -  l)ky  , 
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^41  =  (<t>  +  ^)kz  -  we,^ , 

^42  =  "  u(y  -  l)kz  +  wkx  ,  ^43  =  “  v(y  -  l)kz  +  wky  ,  (B.IO) 

S44  =  (2  -  yjwkz  +  .  S45  =  (y  -  l)kz  , 


Sjj  —  {(4>  +  t|j)  ho}6|^  « 


S52  “  hokx  u(y  1)6^  ,  S53  —  hoky  v(y  1)0)^  , 


(B.ll) 


S54  —  hokz  w(y  l)B|f  t  S55  —  Pic  +  (y  1)®ic  > 

where 

<J)  =  "  2  +  V*  +  w^)  , 

tj)  =  —  h(Y  -  1)  . 


(B.12) 


For  a  perfect  gas,  y  will  be  given  by  the  “standard”  ratio  of  specific  heats  y,  ^  will  be  identically 
zero,  and  the  above  real  gas  formulation  wiU  identically  reduce  to  its  perfect  gas  counterpart 


2.  PRIMITIVE  VARIABLE  JACOBIAN  MATRIX 


In  order  to  work  with  a  more  amiable  eigensystem,  the  conservative  variable  system  is  con- 
voted  to  one  based  on  the  primitive  variables,  where  q  =  J[q,  u,  v,  w,  p]'*^  is  die  vector  of  primi¬ 
tive  variables.  The  transformation  matrix  M,  given  by  M  =  dQ/dq,  is 


M  = 


■  1 
u 

V 

w 

<l> 

Y-  1 

■ 


0  0  0  0 

c  0  0  0 

0  Q  0  0 

0  0  e  0 

CU  ev  QW  ^ 


(B.13) 


Using  die  similarity  transformation,  x  =  M~^?M,the  primitive  variable  Jacobian  matrix  x  is 
given  by 
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‘Pic  Q^X  Oky  ekz  0 

0  Pk  0  0 

*=00  Pi  ®  T  •  ® 

0  -.0  0  Pi  ^ 

0  Ga^kx  Ga%  pa^kj 

which  is  identical  to  the  perfect  gas  formulation  obtained  by  Janus  [Reference  42]. 

The  eigenvalues,  obtained  by  solving  the  equation  be  -  XII  =  0,  are 

Xj[  =  X|  =  X|  =  Pij  , 

X{  =  +  alAkl  ,  (B.15) 

XJ  =  p,,  -  alAkl  . 

Considering  the  equation  (x  -  XI)r  =  0  (where  r  are  the  right  dgenvectors),  the  right  eigenvec¬ 
tor  matrix  P|.  for  the  primitive  variable  Jacobians  is  obtained  as 


okz 

ky 

~  kx 
0 
0 


a  a 
kx  _  ^ 

j2  n 

h  -h 

J2 

kz  _  ^ 

J2  J2 

t?a  a^a 


(B.16) 


where  a  =  Q/&J2  and  the  right  eigenvectors  are  given  by  columns.  Similarly,  the  matrix  of  left 
eigenvectors,  Pj['  ^  can  be  obtained 
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kz 

0 


0 

0 


k  h  k 

J2  J2  J2 

-  ii  _  ^ 

J2  J2  J2 


aa^ 

kz 

aa^  » 


(B.17) 


where  P  =  1  /Qa  J2  and  the  left  eigenvectors  are  given  by  rows.  At  this  point,  the  eigenvector 
matrices  for  the  cmiservative  variable  Jacobians  are  easily  obtained  frmn  the  following  idations 


Tk  =  ^®*k  * 


(B.18) 


where  T,^  is  the  matrix  whose  columns  are  the  right  eigenvectors  and  ^  is  the  matrix  whose 
rows  are  the  1^  eigenvectors.  Perfonning  the  matrix  multiplicaticMis,  the  right  eigenvectors  are 
found  to  be 


ctkx 

aukx 

avkx  +  ekz 
awkx  -  Qky 

oky 

auky  -  Qkz 
ovky 

awky  +  Qkx 

+  C(wkx  -  ukz) 
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r,  = 


ctkz 

auicz  +  Qky 
ovkz  -  ekx 
awkz 


+  o(uky  -  vkx) 


=  a\ 


1 

u  ±  akx 

V  ±  akx 

w  ±  akz 

<|)  +  —  ij) 


Y-  1 


±  a6i 


and  the  left  eigenvectors  are 


li  = 


5(1  -  f 

^ui. 

cta^ 

(T  “  l)..r.  .  kz 


oa-' 


'vkx  +  ^ 


(Y-  Dr 


aa‘‘ 


wkx  -  -j 

(Y-Dr 

aa2  * 


h  ~~ 


<l>  +  ^\r,  ^  1 .  r 


«r  aa2  / 


^\ky  +  i(ukz  -  wkx) 


(Y-  D 


ukv  -  & 


oa-* 


oa-* 


C7-  D 


cca 


2 


(B.19) 


lid 


Similar  to  what  was  seen  for  the  inviscid  Jacabians,  the  eigenvectors  will  nicely  reduce  to  their 
perfect  gas  formulation. 


117 


APPENDIX  C 
SPLUT  FLUX  JACOBIANS 
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The  following  derivation  for  the  split  flux  Jacobians  follows  closely  the  perfect  gas  develop¬ 
ment  given  by  Belk  [Reference  35],  Therefore,  only  the  steps  needed  to  clarify  the  real  gas  ex¬ 
tensions  will  be  given. 

1.  “TRUE”  JACOBIANS 

The  generalized  Steger- Wanning  split  flux  vectors  for  an  arbitrary  gas,  as  developed  by  Vi- 
nolcur  and  Montagnd  [Reference  33],  are 

Q 

Q\1 
QV 

gw 


±  -  r-  1 


Kf=  p 


’(•v.  -  r^r) 


, 


K*  s=  — 

2r 


g(u  +  akx) 
g(v  -I-  aky) 
g(w  +  akz) 

g(ho  +  a0|.) 


If  ±  =  JL 
^5  2r 


g(u  -  akx) 
g(v  -  aky) 
g(w  -  akz) 

Q(ho  -  a0ij) 


where  the  eigenvalues  can  be  written  as 

X  =  JIVkI- 


4  Z 


(C.1) 


The  “true”  (analytical)  Jacobian,  denoted  by  A*,  is  given  by 


3Q  • 


or 


-±  ^  -±  V- 

A  =  I  A.  =  x  ■ 


-l-±  _  -r-± 

i=lv4,5  i*M,5 

The  split  flux  vectors  can  be  redefined  as 

Kt  =  x;*>  gi  K'  , 


(C.2) 


(C.3) 


(C.4) 


(C.5) 
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2.  “X”  CONTRTOUnONS 

The  contributions  to  the  split  flux  Jacobians  arising  from  the  partial  differentiation  of  the 
eigenvalues  can  be  written  as 

^  p  +  ^[1-4  +  Lj] 

u^-^JLx  +  ^[(u  +  akx)L4  +  (u  -  akx)L5] 

+  ^[(v  +  aky)L4  +  (V  -  aky)L5]  (C.12) 

+  ^[(w  +  akz)L4  +  (w  -  ak2)L3] 

ho  —  r  ®®k)^5] 

In  the  above,  the  vectors  L^,  L4  and  L5,  are  the  eigenvalue  derivatives,  which  are  found  to  be 
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(^■L  =  0  • 

("1)22  =  ‘  ’  (“■)23  =  0  -  (“:)24  =  °  •  (“■)25  “  °  ■ 

(C.18) 

0 

II 

»  »-< 

(a;),^  =  0  ,  (a;)^^  =  1  ,  (a;)^  =  0  ,  (a;)^^  =  0  , 

(C.19) 

0 

II 

« 

(>■)«  = « •  (•:)«  = » •  (*;)«  =  ‘  •  (»;)«  =  0  • 

(C.20) 

(*■),.  = 

^P  2ea  aa  .  ea^  ar  a^ 
sq  r  -  idQ  (F  -  1)2 ap  r  -  1  ’ 

=  iL  _  3a  .  dT 

V  V52  aQu  r  -  1  aou  ^  (r  -  i)aaou  ’ 

fa*^  =  iL  _  3a  .  Qa^  gp 

V  V53  agv  r-  laQv  (r-  i)2aQv  * 

fa*^  =  JH-  -  3a  .  Qa^  ar 

V  V54  aQw  r-iapw  (r-i)2aew’ 

U*A  =  JP-  -  ^a  .  Qa^  ar  .  1 

V  V55  agco  r-  lagco  (r-  D^apco 

and  the  components  of  the  matrices,  A44  =  aK4^/aQ,  are 

(^^)ii  =  1’  Mn  “  =  0,  (a:^)^^  =  0,  (a:^)^^  =  0, 


1  ±  ®^*agu  ’  (“^>5)23  -  ±  ekx^  , 

{“^•5)25  "  ’ 


*  ®^*aQW  ’ 


(*♦>3)31  =  ±  aky  ±  Qky||  , 

(*^•5)32  ^  ’  (®^-^)33  ^  ^  ’ 

(^•^)34  "  ®^ya^  ’  ’ 

=  ±  ak.  ±  Qk.||  . 

(®^-^)42  ^  ’  (*^’5)43  ’ 


(C.21) 


(C.22) 


(C.23) 


(C.24) 
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*  *  “  *  ek 


^dQQo  ’ 


(C.25) 


(^’^)51  ~  dQ  ’ 

(^*5)52  ^  ^  ^  ±  akx  , 

(^•5)54  "  ^  ±  akz  . 

(^•5)55  ^  a^  ^^i^aot  ■*■  ^  • 


(C.26) 


5.  SPEED  OF  SOUND  DERIVATIVES 


Recalling  the  definition  for  the  speed  of  sound 

a^  =  r— 

a  i  Q  » 

and  differentiating,  the  following  relation  is  obtained 

M  =  ^iL  +  r  dp  _  ci£ 
aqi  2raqi  Tgadqi  Tgdqi 

Using  Equation  (C.28),  the  speed  of  sound  derivatives  are  found  to  be 

ia  _  a  ar  I  r  a 

dQ  ir  dQ  Tqa  dQ  2q  ’ 

aa  ^  ^JL  +  r 
aQu  2r  aQu  2Qa  agu  ’ 

aa  ^  _a_  jr  .  JLiP_ 
agv  2r  agv  2Qa  agv  ’ 

aa  ^  _a^_ar.  ,  r  ^p 

aQw  2r  a^w  2ea  agw  ’ 

aa  ^  _a__ar.  .  r  ^p 

^QCo  2r  agco  2Qa  ageo  ’ 


(C.27) 


(C.28) 


(C.29) 
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6.  SPLIT  FLUX  JACOBIAN  MATRIX 


The  split  flux  Jacobian  can  be  constructed  by  summing  up  the  contributions  already  pres¬ 
ented.  The  components  of  the  split  flux  Jacobian  matrix  A^  are  given  below,  where  the  follow¬ 
ing  relations  are  used 

a(±) 


Ii  =  ii  . 

(C.30) 

a(±)  a(±) 

I2  -  I4  “  l5 

(C.31) 

a(±)  a(±) 

I3  —  I4  +  Is  A 

(C.32) 

a(±)  j  /a(±) 

14  =  ii  -5(i4 

a(±)\ 

•f  is  1  . 

(C.33) 

The  components  of  the  first  row  are 


a  =  icl  —  ^I  -1-  -^I  +  C  —  -I-  C  — 
^11  ^  ^  ^  ’ 


^12  - 


dQ 

+  c 


Eia  —  kwli  Cl 


dp 


+  c 


iBapu  ’ 

ar 


‘13  -  p  *4  -  ^  ''IB  aQv  ’ 


ai4  -  kzli  i^^4  +  Cia^ 


*®3qw  ’ 


a  =  k  I  —  -I-  -iLr  +  c  —  -t-  C 

“11  *^t^l  pi4  T  ^pl2  ^  '-lA  A*  T 


^12  “  ^x^l  p  I4  +  ^IB 


lAag  ^  ^IB^q  » 

ar 

agu  ’ 


*13  ~  ^y^l  'F^4  ^lA 


aQu 

iL  +  c 

agv  aQv  ’ 


^zt  j_  /->  ar 


‘y*!  p 

5i4  =  kJi  -  -^14  +  CiA^^  +  CjB 
*15  “  C 

where 


JP_  +C 

'lAa^eo  ^i^ageo  ’ 


a^w  ’ 


c  =  -^i 

'-lA  4^12  > 


(C.34) 


(C.35) 


124 


c  =  — 

'-IB  -  p2 


''IB  ^  p2[M4  4^2]  • 

The  components  of  the  second  row  are 

=  -  uS^Ij  -  ^(u  +  kxP,^)l2  +  ^kxktl2  +  +  ^2b|q  ’ 

Ml  -  — K-'U  +  #(1  +  kx')lj  + 


(C.36) 


=  (>*x  +  M,  -  — +  A(1  + 

iilr  <alr  Ir 

523  =  ^yh  -  y  1 

124  “ 


+  C 

^  ^2Ba0u  . 


=  C2A^^  +  QbI^  . 
^^dgco  “3Qeo 


®25  - 
where 


+r  ^  +C 

2p  I2  +  ^2A3Qy  +  ^2BSQy  • 

^xkzr  1  p  I  p  3r 

2r  ^2  ^aQw  • 

ar 

ape©  ’ 

(C.37) 

~  ^[  (tt  +  kxP)j)l3  +  2akxl2]  » 

(C.38) 

^xPk)^2j  • 

(C.39) 

Ul  UlC  UIUU  lUW  019 

=  -  V0,jli  jr^l4  -  ^(v  +  kyp,j)l2  +  ^kyktl2  +  +  Cgel^  , 

=  vie  I  —  +  C,  +  c, 

vKxii  p  I4  +  2r  ^  ^^dQu  ^  ^^dgu  ’ 

I„  -  (vk,  +  Px)I.  -  +  A(1  +  k,\  +  . 

^  ..C  a.  .n 


»  SS  C  4.  C 

533  +^3Bdgeo  ' 

where 


531 

532 

*33 

534 


+  +C,B/r'  . 

2r  ^  ^^aQw  ^®aQw 

ar 

^Q®o 

(C.40) 

^3A  ”  +  ^yPk)^3  ■•■  2akyl2]  , 

(C.41) 

^3B  “  ”  4(^  ■*■  ^yPk)^2]  • 

(C.42) 
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The  components  of  the  fourth  row  are 


5  -  wt  T  4.  akxkz.  .  ^  ^p  4.  r  5F 

a42  wkxlj  p  I4  +  2r  ^2  +  Qa^qu  ■*'^4BaQu  ’ 

wkv  akvkr  do  ;5r 

143  wkyii  p  I4+  2r  ^2  +  C4Aagv  ■‘■^«aev  ’ 

~  wk,  +  6v  a  ~  2 

S44  =  (wkj  +  Pt)Ii - V-“U  +  ^(1  +  k.  )Ij  +  C,^ 

ap  .  f-  sr 

dQW  4B  aQw  ’ 

a  -  C  +  C 

^45  ’-4Aageo  ^*aoeo  ’ 

(C.43) 

where 

^4A  “  ■*■  k2Pjt)l3  +  2akxl2]  , 

(C.44) 

^4B  =  ^  [wpkl4  -  f  (w  +  ^A)h]  • 

(C.45) 

The  components  of  the  fifth  row  are 


~  ~  2 
*51  =  -  Mfcii  +  *  *1«  -  ^(h,  +  e,p,)ij  -  ^k,k,ij 

—  f~'  C' 

'5AaQ  +  '-sBag 


a  “  h  k  T  —  T  4.  4>  ft  M  —  C  ^P  4.  c 

*52  “  hoKxli  p  ^  2r  ^  ^5AaQu  ^  ^5Ba^  » 

I5,  =  h4,I.  -  +  ^(9^  +  WIj  -  C,*^  +  C„^  , 

*«  =  hjij,  -  +  ^(0^  +  p^)Ij  -  +  Cb^  . 

*35  =  Wl  “  ^14  +  “  '-SAj^  +  Cto-^ 


agco  ’ 


whm« 


CjA  ~  1^14  "^^3  ■*■  '*'  ®kPk)^2  ■•■  ^^2  > 

^5B  ~  j^^^Pk^4  "*■  2^  ^  ®k^k)^2j  * 


(C.46) 


(C.47) 

(C.48) 
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The  pressure  derivatives  appearing  in  the  above  formulation  are  given  by  Equations  (A.22)  - 
(A.26) .  In  the  present  work,  the  “gamma”  derivatives  are  assumed  to  be  zero,  since  their  deter¬ 
mination  is  much  involved  and  would  necessitate  the  use  of  numerical  Jacobians,  either  fully  or 
in  some  form  of  hybrid  scheme.  The  above  formulation  differs  from  the  perfect  gas  formulation 
in  the  evaluation  of  the  pressure  derivatives,  tire  use  of  the  isentrc^ic  index  Fas  opposed  to  the 
ratio  of  specific  heats  y,  and  the  q)pearance  of  the  “gamma”  derivatives.  With  these  consider¬ 
ations  in  mind,  the  split  flux  Jacobian  matrix  for  a  gas  in  chemical  equilibrium  can  be  shown  to 
reduce  to  its  prefect  gas  counterpart. 
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APPENDIX  D 

ROE  AVERAGED  VARIABLES 
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1.  JUMP  CONDITIONS 


The  key  step  in  the  development  of  an  approximate  Riemann  solver  is  to  determine  averages 
of  eigenvalues,  eigenvectors,  and  wave  strengths  such  that  the  jump  conditions 


(D.l) 


ISI  =  ^6,  Xi  fi , 


(P-2) 


are  satisfied.  The  averaged  eigenvalues  read 


AAA 


Xi  =  Xj  =  X,3  =  JlVkiPjj  , 

X4  =  JlVkK  Pk  +  a  )  . 


(D.3) 


X3  =  JIVkl(  Pk  -  a  )  . 

where  the  directional  subscript  k  has  been  drqpped  for  clarity.  The  averaged  eigenvectors  read 


avkx  +  Qkz 
awkx  -  pky 

-  a^/(Y  “  1))  +  C(vks  -  wky) 


A  A«  A« 

auky  -  Qkz 


awky  +  Qkx 

aky|ho  -  a^/(Y  -  1))  +  e(wkx  “  ukx) 


u±ak. 


=  a  V  ±  akx 
w  ±  akz 

A  A  A 

ho  ±  aOk 


(D.4) 
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and  the  wave  strengths  read 


A  ky  j 

[m-M] 
V  a  ) 

1  +  kzlv]]  -  kylw]]  , 

Q> 

10 

II 

\  a  y 

1  -  yul  +  kxlwl  , 

[m-M] 

^  a  J 

1  +  kyM  -  kzlvll  , 

II 

'  ^  ±  le  J 

{  Qti 

) 

2.  IQI- EQUATIONS 

The  system  givoi  by  Equation  (D.l)  can  be  writtra  as 

lel  =  ffeJ . 

Icull  =  uIqI  +  ^lul  , 

IqvII  =  vttej  +  civi , 

IqwII  =  wIqI  +  qIIwI  , 

jk  ^ 

ttOCoB  =  — —IqI  +  7-^— -llpll  +  ^u|[u]l  +  vlvj  +  wlwl)  . 

Y  —  1  Y  ~  1 

3.  |[FI|  -  EQUATIONS 

The  system  given  by  Equatitxi  (D.2)  can  be  written  as 

ttePkl  =  Pklloll  +  My . 

loujk + py  =  upijUei + cPkM + Miy + kxipii , 


(DJ) 


(D.6) 

(D.7) 

(D.8) 

(D.9) 

(D.IO) 


(D.ll) 

(D.12) 
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(D.14) 


tlQw3k  +  pkj]  =  wPj^iieii  +  ePkII'^1  +  Qwie j  +  kjp]] , 

*  A 

lehoPk  -  pkj  =  A  _  ^  PkfeU  +  0Pk(“II“ll  +  vlvl  +  wlwl) 

A  A  2  A 

+  T^W  +  e,llp]l  +  r'^iej  . 

Y  -  1  Y  “  1 


(D.15) 


Multiplying  each  IQl-momcntum  equation.  Equations  (D.7)  -  (D.9),  by  their  correspond¬ 
ing  directional  metric,  as  follows 


kxlloull  *  ukxlloll  +  e^xllull  . 

(D.16) 

kylQvIl  =  Vkylel  +  Qkyllvll  , 

(D.17) 

kxiowll  =  wkxlell  +  Qkxllwll  , 

(D.18) 

and  adding  all  diree  equations,  the  following  is  obtained 

IkxQuI  +  IlkyQv])  +  Ikxewll  =  (ukx  +  vky  +  wki)|lQ]l 

+  e(|Iukx  +  vky  -I-  wkj)  . 


Adding  [QHlq  to  both  sides  and  regrotqnng  yields 

IcPkl  =  Mel  +  MM  *  (D^O) 

which  is  identically  Equation  (D.ll).  Following  the  same  procedure  for  the  |[F]]-m(Hnentum 
equatimis.  Equations  (D.12)  -  (D.14),  and  adding  |[QP]^]]kt  to  both  sides,  die  result  is 


ttcPkl  =  Pkttell  +  2cp^|iej  .  (D.21) 

However,  noting  diat  the  jump  in*  the  temporal  metric  is  zero,  the  jump  in  the  relative  contravari- 
ant  velocity  can  be  written  as 


IPkll  =  ttukx  +  vky  +  wk J  -I-  IkJI  =  10  J  .  (J>22) 

Now  Equadcms  (D.20)  and  (D.21)  can  be  rewritten  as 


IcPkll  “  Pklell  +  MPkl  * 


(D.23) 


(D.24) 


IqPicI]  ~  Pkttol  2qPj^IIP  j  . 

A  * 

Solving  for  q  and  the  following  averages  are  obtained 

Q  =  =  JqTqi  ,  (D.25) 

A  (Pk^)  _  (Pk)r>/^  +  (Pk)l'/^  /n'rr\ 

'’‘■W  ^  +  ^  ' 


where  the  arithmetic  average,  <  *>=[(*  )r  +  ( ’  )i]/2,  has  been  utilized.  Similar  to  die  result 
for  the  relative  ccmtravahant  velocity,  the  averages  for  the  absolute  contravariant  velocity  and 
Cartesian  velocities  are  given  by 


0,j  =  31»(0k)  ’  “  ”  ^  ^  ~  9l(w)  , 

(D.27) 

where  denotes  the  “standard”  Roe  average 

aaa  X  _  (’)T'lQt  +  i')\y/Qi 

(D.28) 

A 

Looking  at  the  energy  equations.  Equations  (D.IO)  and  (D.IS),  multiplying  the  former  by 
and  subtracting  the  result  from  the  latter,  die  following  is  obtained 

A  A  +  tb 

llQMk  “  pM  PkIIc«oI  =  0kllp]l  +  A  _  j  '^lOkl  • 

(D.29) 

Recalling  that 

A  a2  a 

tjj  =  a  -  (Y  -  l)h  , 

(D.30) 

Jk  Y  -  1/a2  ^  a2  ^  a2\ 

<|)  =  “f"!"  ^  ^  j  ’ 

(D.31) 

iPk^  ”  » 

(D.32) 

Equation  (D.29)  can  be  rewritten  as 

leMkl  =  ehottPkl  +  PklloM  .  (D33) 

A 

Solving  the  above  results  in  the  following  average  for  the  total  enthalpy  ho 
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ho  =  ^R>(ho)  .  (D.34) 

Returning  to  Equation  (D.IO),  and  splitting  the  total  energy  into  internal  and  kinetic  contribu¬ 
tions,  results  in  the  following 

A  A 

loeol  =  I  _  "Hell  +  ,  (D.35) 

From  Equation  (D.24)  the  following  formulas  for  the  Cartesian  velocities  can  be  obtained 

Hqu^J  =  u^Io]]  +  2qu|Iu]]  , 
lev^I  =  +  2ev|v]l  , 

lew^I  =  w^|[q]1  -I-  2qw|[w]1  . 


Substituting  the  above  relations  into  Equation  (D.35),  the  following  formulation  for  the  internal 
energy  is  obtained 

A 

M  =  -  T-^lIell  .  (D.36) 

Y  -  1  Y  ~  1 

In  order  to  obtain  all  jumps  in  terms  of  q  and  T,  the  following  relations  are  introduced 


IpI  =  eRlTll  +  RTIqI  +  QT^^Rill^]]  . 

i=l 

leel  =  e  idjlej  +  e  +  i(|)e|ell  • 

i=l'  '  i=l  i=l'  ' 


(D.37) 

(D.38) 


In  the  above,  averages  for  mixture  gas  constant,  temperature,  species  internal  energy  and  species 
mass  fraction  are  given  by  the  “standard”  Roe  formulas,  as  follows 


R  =  3l(R)  ,  t  =  3l(D  ,  Ci  =  ,  Yj  =  ^(Y^)  . 

The  jump  in  species  internal  energy  can  be  written  as 

lej  =  CvjJT]]  , 


where  is  an  integral  average 


Cv.(x)dt  . 


(D.39) 

(D.40) 


(D.41) 
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The  jump  in  the  species  mass  fractions  can  be  obtained  as 


l|l  =  I  (^)de  +  |'(§)dT  ,  (D.42) 

and  rewritten  using  integral  averages  as  follows 

=  (W)m  +  (^)lTll  .  (D.43) 

Vinokur  [Reference  7]  developed  a  procedure  that  would  identically  satisfy  Equation  (D.43). 
However,  this  procedure  is  complicated  and  numerically  expensive.  The  simplest  approximation 
for  the  integral  averages  is  the  trapezoidal  rule,  i.e. 


(D.44) 


which  is  the  choice  adopted  in  this  woilc.  Substituting  Equations  (D.37)  -  (D.43)  into  Equation 
(D.36)  and  combining  terms,  results  in  the  following 


=  0  .  (D.45) 


In  order  for  Equation  (D.45)  to  be  satisfied,  both  coefficients,  denoted  by  {  },  must  be  identical¬ 
ly  zero,  thus  giving  tiie  following  two  equations 


(D.46) 

(D.47) 


Solving  Equation  (DM)  results  in  the  average  for  the  ratio  of  enthalpy  and  internal  energy  de¬ 
rivatives,  as  follows 
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Y  =  1  + 


i=l 


(D.48) 


where 


Cv  “ 


(D.49) 


Equation  (D.47)  can  be  solved  for  tp,  as  follows 


A  ^  A  A 

\|)  =  FRT  +  (y 


(D.50) 


where  the  average  for  the  isentrq)ic  index  is  taken  to  have  the  same  fonn  as  given  in  Equation 
(A.  16),  as  follows 


(D.51) 


Substituting  Equation  (D  JO)  into  Equation  (D.30)  results  in  the  following  average  for  the  speed 
of  sound 


A  2  a  fL 

a  =  FRT  +  (Y  -  1) 


/»2  N 

A  q 


A  n  A  ,  A  A 

hc-T-ZY,4i-RT 


i=l 


(DJ2) 


which  can  be  rearranged  in  the  form 

a21 


a2  a  -v 

a  =  (y  -  1) 


A  Q 

ho  ~  ^1  + 


|(r  -  Y  +  I)RT  -  -  »  Z  , 


(DJ3) 


2  2  2  2  A 

whereq  =u  +  v  +  w  .  For  a  perfect  gas,  where  T  =  y  =  Y.  the  averages  for  the  Cartesian 
velocity  components  and  total  enthalpy  are  the  same,  and 


=  (D^) 

i-l  ' 

Consequently,  the  term  in  braces  in  Equation  (D  J3)  will  be  identically  zero,  and  the  relation  for 
the  averaged  speed  of  sound  will  reduce  to 


a2 

a 


(Y-  1) 


which  is  the  original  formula  developed  by  Roe  [Reference  34]. 


(D.55) 
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APPENDIX  E 

DISaiETIZED  VISCOUS  JACOBIANS 


1 .  DISCRETIZATION  OF  VISCOUS  JACOBIANS 

The  following  are  the  components  of  the  discretized  viscous  Jacobian  matrix,  taken  at  cell 
face  1  +  1/2 

i^y)n  =  ° 

=  =F  lAkljj^ll  +  |k,  . 

(.*)„=  ^IAld{[ttk.k,]_^|>)J, 


(5.^)3,  -T  IAkl|[Hk.k,]^^|i)^J  , 

(5,*)33  =  T  IAkl|[M  +  . 


(E.2) 
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(E.3) 


(E.4) 
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+  ^IAkl 


Jl  +  : 


(E.5) 


In  the  above. 


V  ^xaic  +  *^yak^^"ak  ’ 


(E.6) 


q2  =  , 


(E.7) 


and 


(*)i+  —  (‘)i  »  (’)i-  —  (')i+i  •  (E.8) 

The  temperature  derivatives  can  be  derived  in  a  similar  fashion  as  the  pressure  derivatives,  see 
Appendix  A.  The  final  results  read 
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